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Abstract 

This  work  presents  a  novel  approach  to  code  phase  multipath  mitigation  for  Global  Posi¬ 
tioning  System  (GPS)  receivers.  It  uses  the  power  and  complex  cepstra  for  multipath  detection 
and  mitigation  prior  to  code  phase  tracking  by  a  standard  non-coherent  delay  lock  loop.  Cepstral 
theory  is  presented  to  demonstrate  how  multipath  reflection  delays  can  be  detected  through  the  use 
of  the  power  cepstrum.  Filtering  can  then  be  performed  on  the  complex  cepstrum  to  remove  mul¬ 
tipath  effects  in  the  cepstral  domain.  Finally,  an  inverse  complex  cepstrum  is  calculated  yielding  a 
theoretically  multipath  free  direct  path  estimate  in  the  time  domain.  Simulations  are  presented  to 
verify  the  applicability  of  cepstral  techniques  to  the  problem  of  GPS  multipath  mitigation.  Results 
show  that,  under  noiseless  conditions,  cepstral  processing  prior  to  code  tracking  by  a  standard  non¬ 
coherent  delay  lock  loop  leads  to  lower  code  tracking  biases  than  direct  tracking  of  the  composite 
multipath  signal  by  a  narrow  correlator  receiver.  An  exception  to  this  general  rule  occurs  at  a  mul¬ 
tipath  reflection  delay  of  exactly  1.0  chip,  relative  to  the  direct  path,  where  the  cepstral  processing 
provides  no  improvement.  Additionally,  cepstral  processing  provides  little  or  no  improvement  over 
a  narrow  correlator  non-coherent  delay  lock  loop  when  the  multipath  delay  is  small,  approximately 
0.1  chip  or  less.  However,  this  weakness  is  common  for  other  multipath  mitigation  techniques  as 
well.  Finally,  this  work  shows  that  cepstral  processing  is  highly  sensitive  to  additive  white  Gaussian 
noise  effects,  leading  to  the  conclusion  that  methods  of  limiting  noise  effects  must  be  developed 
before  this  technique  will  be  applicable  in  actual  GPS  receivers. 
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Cepstral  Processing  For  GPS  Multipath  Detection  and  Mitigation 

L  Introduction 

1.1  Overview 

Satellite  navigation  systems  have  ushered  in  a  new  era  in  precision  navigation.  The  Global 
Positioning  System  (GPS)  in  particular  provides  highly  accurate  position  solutions  for  both  civilian 
and  military  applications.  This  has  led  to  it’s  use  in  every  arena  from  personal  recreation  to  highly 
precise  “smart”  bombs.  However,  a  primary  error  source  still  exists  in  GPS,  ranging  errors  caused 
by  multipath  interference.  In  fact,  multipath  signals  are  the  dominant  error  source  for  Differential 
GPS  (DGPS).  Thus  many  efforts  focus  on  developing  techniques  to  mitigate  multipath  errors  in 
DGPS. 

1.2  GPS  Background 

The  Global  Positioning  System  is  a  direct  sequence  spread  spectrum  (DS/SS)  satellite  nav¬ 
igation  system  consisting  of  24  satellites  providing  navigation  information  to  passive  receivers. 
The  system  operates  using  the  principle  of  trilate ration.  Within  the  received  signal,  the  satellite 
transmitter  embeds  precise  timing  information.  The  receiver  uses  the  signal  propagation  delay  to 
calculate  the  distance  to  the  satellite  using  d  =  Vp  ^  {U  -  U)  where  d  is  the  distance  from  the 
receiver  to  the  satellite,  Vp  is  the  velocity  of  propagation  (usually  assumed  to  be  the  speed  of  light 
in  a  vacuum),  tt  is  the  time  of  transmission,  and  tr  is  the  time  of  receipt. 

In  addition  to  transmitting  precise  timing  information,  each  satellite  transmits  ephemeris  data 
which  provides  the  satellite’s  position.  From  knowledge  of  the  positions  of  at  least  four  satellites 
and  the  range  to  each  satellite,  a  user  position  solution  is  calculated.  This  still  leaves  two  important 
questions  unanswered.  First,  how  does  the  receiver  distinguish  between  satellites?  Second,  how 
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does  the  receiver  determine  the  precise  time  of  transit  for  the  signal  being  transmitted  from  the 
satellite?  Both  of  these  questions  will  be  covered  in  the  sections  that  follow. 

1.2.1  The  GPS  Signal  Before  exatnining  the  GPS  signal,  consider  a  general  direct 
sequence  spread  spectrum  signal.  This  signal  is  represented  as 

s(t)  =  cos{uJot  Oc{t)  +  Od{t))  (1) 

where  s{t)  is  the  transmitted  signal,  P  is  the  carrier  power,  loq  is  the  carrier  angular  frequency,  Oc 
is  the  spreading  code,  6d  is  the  data  stream.  Assuming  antipodal  binary  phase  shift  keying  (BPSK) 
is  used  for  the  spreading  code  and  data  modulation,  an  equivalent  form  of  Equation  1  is  given  in 
Equation  2 

s{t)  =  V^c{t)d{t)  cos{u)ot)  (2) 

where  d{t)  is  the  data  stream  and  c{t)  is  the  spreading  code. 

The  Global  Positioning  System  consists  of  satellites  transmitting  two  distinct  DS/SS  signals. 
The  first  signal  is  called  Precision  code  (P-code).  This  signal  is  encrypted,  a  conversion  to  Y- 
Code,  before  transmission,  and  is  reserved  for  military  users.  The  P-code  signal  provides  position 
estimates  with  accuracy  on  the  order  of  10  m  in  the  absence  of  multipath  interference.  The  second, 
Coarse/Acquisition  (C/A)  code,  is  transmitted  for  all  users.  The  C/A  code  signal  is  subjected  to 
selective  availability  (SA),  a  form  of  intentional  signal  degradation.  With  SA  in  operation,  C/A 
code  provides  position  estimates  with  accuracy  on  the  order  of  100  m  in  the  absence  of  multipath 
interference.  P-code  is  transmitted  at  two  carrier  frequencies,  LI  and  L2;  whereas  C/A  code  is 
available  only  on  LI.  The  properties  of  the  GPS  signal  are  summarized  in  Table  1  (5). 

1.2.2  Spreading  Code.  The  Global  Positioning  System  uses  orthogonal  Gold  codes  to 
spread  the  signal  spectrum.  For  a  discussion  of  Gold  codes  see  (9,  10).  These  codes  allow  GPS 
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Table  1  GPS  Spreading  Code  Characteristics 


Parameter 

C/A  Code 

P  Code 

Data  Rate 

50  Hz 

50  Hz 

Chip  Rate  =  1/Tc 

1.023  Mchips/sec 

10.23  Mchips/sec 

Code  Period  =  N  (chips) 

1023  Chips  (1  ms) 

6  X  10^^  Chips  (1  week) 

Carrier  Band  Designation 

Ll 

1575.42  MHz 

Ll,  L2 

1575.42,  1227.6  MHz 
respectively 

to  utilize  code  division  multiple  access  (CDMA)  to  distinguish  which  satellite  transmitted  a  given 
signal.  Due  to  the  orthogonality  of  the  codes,  the  correlation  of  the  received  signal  will  be  essentially 
zero  if  the  locally  generated  code  is  not  the  same  as  that  of  the  received  signal. 

GPS  receivers  exploit  knowledge  of  the  spreading  code  correlation  function  to  synchronize  a 
locally  generated  replica  of  the  spreading  code  to  the  received  spreading  code.  This  synchronization 
is  needed  to  determine  range  to  the  satellites;  additionally,  this  tracking  of  the  spreading  code 
accomplishes  despreading  of  the  received  signal.  The  code  autocorrelation  function,  Rc{t)^  is 
defined  as 

1 

^  'nTc  Jo  4*)4t+TTc)dt 

where  r  is  the  independent  (time  shift)  variable,  NTc  is  the  code  period  in  seconds,  and  c(-)  is  the 
spreading  code. 

For  maximal  length  pseudorandom  noise  (PRN)  codes,  the  fundamental  period  of  the  code 
autocorrelation  function  is  (5,  11) 


Rc{r)  =  < 


+  kl  <  1 

-i  l<|r|<(iV-l) 

:|r  |-(7V-l)](l  +  i)-i  (iV-l)<|r|<iV 


(4) 


where  N  is  the  code  period  in  chips. 
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Since  iV  1,  Equation  4  can  be  approximated  as 


Rc{r)  «  { 


1-  |t| 
0 


kl  <  1 

elsfAvhere 


(5) 


A  plot  of  the  Equation  5  is  shown  in  Figure  1.  Although  Figure  1  shows  only  a  single  peak, 
the  periodicity  of  the  spreading  code  results  in  a  periodic  autocorrelation  function.  Therefore,  the 
peak  of  Figure  1  is  repeated  every  N  chips. 


Figure  1  Maximal  Length  Spreading  Code  Autocorrelation  Function 


L2.3  Code  Tracking.  As  alluded  to  earlier,  code  tracking  in  a  DS/SS  receiver  is  accom¬ 
plished  by  computing  the  correlation  of  the  received  and  locally  generated  spreading  codes.  The 
locally  generated  code  is  shifted  in  time  until  the  correlation  is  maximized.  Maximization  of  the 
correlation  function  occurs  at  r  =  0.  Referring  to  Equation  3,  it  is  seen  that  r  =  0  implies  that  the 
locally  generated  code  and  the  received  code  are  synchronized.  By  measuring  the  shift  of  the  local 
code  necessary  to  maximize  the  correlation  function,  the  receiver  can  estimate  the  signal  propaga- 
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tion  delay.  As  discussed  earlier,  the  propagation  delay  is  then  used  to  estimate  the  range  to  the 
satellite. 

l.S  Multipath  Interference 

A  multipath  signal  is  a  reflection  of  the  direct  path  signal.  The  multipath  phenonmenon  is 
illustrated  in  Figure  2. 


Figure  2  The  Multipath  Signal 


From  Figure  2  it  is  seen  that  the  multipath  signal  traverses  a  greater  distance  than  the  direct 
path  signal;  hence,  the  multipath  signal  is  a  delayed  version  of  the  direct  path  signal.  Because 
it  is  assumed  that  the  multipath  signal  is  also  attenuated  with  respect  to  the  direct  path  signal, 
Equation  2  can  be  used  to  represent  the  multipath  component  as 

7ni{t)  =  ai  s{t  -  ~  ai  c{t  -  [2c^o  {t  -  (6) 
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where  ai  is  the  multipath  attenuation  coefficient  and  Trm  is  the  multipath  delay  relative  to  the 
direct  path  signal. 

Using  the  superposition  of  Equations  2  and  6,  the  total  received  signal  is  represented  by 
Equation  7 

M 

r{t)  =  s{t)  +  m{t)  =  c{t)d{t)  cos  (o^o  0  ~  -  Tmi)  \/^  cos[2a;o  {t  -  r^i)] 

i=\ 

(7) 

where  the  summation  indicates  M  >  1  total  multipath  reflections. 

l.S.l  Code  Tracking  in  the  Presence  of  Multipath.  When  code  tracking  in  the  presence 
of  multipath  is  performed  in  the  same  manner  as  code  tracking  in  the  absence  of  multipath,  the 
multipath  term(s)  in  Equation  7  cause  a  distortion  of  the  autocorrelation  function.  This  degrada¬ 
tion,  which  will  be  discussed  later,  leads  to  a  non-zero  tracking  bias  with  respect  to  the  direct  path 
signal.  In  other  words,  rather  than  tracking  the  direct  path  signal  perfectly,  the  multipath  compo¬ 
nents  cause  some  tracking  error  in  the  receiver.  These  effects  will  be  discussed  more  thoroughly  in 
Chapter  III. 

1.4  Previous  Efforts  in  GPS  Multipath  Mitigation 

Efforts  to  eliminate  multipath  errors  in  GPS  receivers  have  focused  on  three  main  areas: 
pre-receiver,  receiver-internal,  and  post-receiver.  Pre-receiver  techniques  typically  include  antenna 
designs  that  minimize  the  antenna  gain  in  the  expected  direction  of  the  multipath  signal.  Alter¬ 
natively,  the  antenna  can  be  designed  to  have  a  physical  barrier,  such  as  a  ground  plane,  blocking 
the  multipath  signal.  Such  antenna  designs  work  well  in  static  applications  where  the  trajectory 
of  the  multipath  signal  is  relatively  constant  and  can  be  estimated.  However,  they  do  not  perform 
well  in  dynamic  situations  where  the  multipath  environment  is  changing  rapidly.  Receiver-internal 
techniques  typically  consist  of  digital  signal  processing  techniques  designed  to  eliminate  the  mul- 
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tipath  components  from  the  received  signal,  or  to  minimize  the  negative  effects  of  the  multipath 
components.  These  techniques  perform  well  even  in  highly  dynamic  multipath  environments.  This 
thesis  fits  into  the  receiver-internal  category.  A  brief  summary  of  other  receiver-internal  design  and 
signal  processing  techniques  previously  developed  follows. 

1.4.1  Narrow  Correlator  Spacing.  The  use  of  a  non-coherent  delay  lock  loop  (NCDLL) 
with  a  more  narrow  correlator  spacing  was  one  of  the  first  receiver  designs  for  multipath  mitigation 
(12).  This  technique  employs  a  code  tracking  delay  lock  loop  with  correlators  spaced  0.1  of  a 
chip  duration  apart  in  time  rather  than  the  standard  1  chip  width.  The  narrow  spacing  results  in 
the  delay  lock  loop  tracking  a  more  narrow  portion  of  the  correlation  peak.  This  leads  to  better 
performance  in  a  multipath  environment  with  majclmum  tracking  errors  being  reduced  by  a  factor 
of  approximately  10  with  respect  to  the  standard  correlator  structure.  This  technique  is  currently 
used  in  many  production  GPS  receivers  for  standard  positioning  system  (SPS)  (i.  e.  C/A  code)  use. 

1.4-2  Multipath  Estimating  Delay  Lock  Loop.  The  multipath  estimating  delay  lock  loop 
(MEDLL)  provides  improvements  over  the  narrow  correlator  spacing  receiver  (13).  In  the  MEDLL, 
the  multipath  signal  amplitudes  and  delays  are  estimated,  and  the  multipath  signal  is  removed 
through  digital  signal  processing.  The  MEDLL  is  implemented  using  a  bank  of  correlators,  each 
delayed  relative  to  the  other  by  a  fraction  of  a  chip.  By  correlating  the  received  signal  using  this 
bank,  a  sampled  version  of  the  distorted  correlation  function  is  generated.  The  MEDLL  then 
uses  a  maximum  likelihood  estimator  to  estimate  the  multipath  signal  parameters.  Finally,  using 
digital  signal  processing,  the  estimated  multipath  components  are  removed  from  the  received  signal. 
This  proprietary  technique  provides  significant  improvement  over  narrow  correlator  spacing  and  is 
implemented  in  some  production  GPS  receivers. 

1.4.3  Modified  RAKE  Delay  Lock  Loop.  The  modified  RAKE  delay  lock  loop  (MRDLL) 
was  proposed  in  an  Air  Force  Institute  of  Technology  thesis  by  Mark  Laxton  (5).  The  MRDLL  has 
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similarity  to  the  MEDLL,  using  a  bank  of  correlators  and  a  maximum  likelihood  estimation  unit 
to  estimate  the  multipath  parameters.  The  MRDLL  employs  an  adaptive  loop  controller  to  adjust 
the  loop  filter  gain,  maintaining  a  fixed  linearized  loop  natural  frequency  and  damping  ratio.  This 
technique  has  a  non-zero  tracking  error  in  the  absence  of  multipath  interference. 

Correlator  Reference  Waveform  Design.  Correlator  reference  waveform  design  is 
one  of  the  latest  techniques  to  be  suggested  for  GPS  multipath  mitigation  (14).  This  technique 
recommends  correlating  the  received  signal  with  a  waveform  designed  explicitly  to  reject  multipath 
signals.  Rather  than  using  a  replica  of  the  GPS  spreading  code  for  correlation  in  the  delay  lock  loop, 
Weill  recommends  using  the  second  or  fourth  derivative  of  the  spreading  waveform.  This  technique 
claims  to  reduce  multipath  ranging  errors  by  increasing  the  range  resolution  of  the  correlation 
process  through  the  use  of  novel  reference  waveforms.  The  advantage  of  the  second  and  fourth 
derivative  correlators  is  that  they  can  be  implemented  using  two  correlators,  one  for  code  tracking, 
and  one  for  acquisition,  data  recovery,  and  data  removal.  This  is  opposed  to  a  bank  of  many 
correlators  for  the  MEDLL  and  the  MRDLL.  Additionally,  this  technique  can  be  implemented 
easily  without  significant  additional  signal  processing,  such  as  that  required  by  the  MEDLL  and 
the  MRDLL.  The  disadvantages  are  that  this  technique  requires  a  slightly  higher  input  signal  to 
noise  ratio  (SNR)  than  the  other  techniques,  and  this  technique  cannot  completely  remove  the 
multipath  interference.  Rather,  correlator  reference  waveform  design  greatly  reduces  the  error 
caused  by  the  multipath  signals.  Like  the  MRDLL,  this  technique  is  very  new,  and  not  likely  to 
have  been  implemented  in  a  production  model  GPS  receiver. 

1.5  Problem  Statement 

The  complex  cepstrum  is  a  tool  originally  investigated  for  echo  determination  in  seismic 
applications.  By  virtue  of  the  fact  that  the  cepstrum  is  a  complex  quantity,  signal  magnitude  and 
phase  information  is  preserved  through  the  application  of  the  complex  cepstrum.  The  signal  can 
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subsequently  be  reconstructed  through  an  inverse  cepstrum  process.  Since  multipath  signals  are 
“echoes”  of  the  direct  path  signal  (see  Equation  6),  the  complex  cepstrum  will  be  investigated 
as  a  tool  to  remove  the  multipath  interference  in  a  composite  signal.  This  thesis  will  pursue  the 
complex  cepstrum  as  a  means  to  sepai’ate  the  multipath  interference  from  the  composite  GPS 
signal.  The  multipath  interference  will  then  be  removed  by  filtering  (often  called  littering)  in  the 
cepstral  domain.  Finally,  the  signal  will  be  transformed  back  into  the  time  domain  for  tracking 
by  a  standard  GPS  receiver.  With  the  multipath  interference  thus  removed,  the  induced  tracking 
error  will  be  removed,  allowing  accurate  tracking  of  the  direct  path  GPS  signal. 

1.6  Objectives 

This  thesis  proposes  a  GPS  receiver  employing  a  complex  cepstrum  process  for  multipath 
removal,  followed  by  a  standard  delay  lock  loop  for  code  tracking.  This  new  design  modifies 
current  GPS  receivers  by  employing  front  end  digital  signal  processing  for  the  removal  of  multipath 
interference.  The  rest  of  the  design  is  a  standard  GPS  receiver. 

The  objectives  of  this  thesis  are: 

1.  Develop  a  complex  cepstrum  filtering  technique  for  GPS  multipath  interference  removal. 

2.  Modify  a  standard  delay  lock  loop  by  adding  the  filtering  to  the  loop  input. 

3.  Characterize  the  multipath  induced  tracking  error  for  the  modified  delay  lock  loop  in  a 
noiseless  environment. 

4.  Characterize  the  multipath  induced  tracking  error  for  the  modified  delay  lock  loop 
operating  with  typical  GPS  signal  to  noise  ratios. 

5.  Compare  and  contrast  the  complex  cepstrum  receiver  performance  to  that  of  a  narrow 
correlator  receiver. 
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1,7  Assumptions 


For  this  thesis,  the  following  assumptions  are  made: 

1.  The  received  signal  consists  of  a  direct  path  and  a  single  reflected  signal. 

2.  The  receiver  has  completed  the  signal  acquisition  phase  prior  to  multipath  interference 
occurring. 

3.  The  GPS  signal  has  been  down  converted  to  an  intermediate  frequency. 

4.  The  received  signal  is  due  to  only  one  GPS  satellite  (i.e.  one  signal  per  channel). 

5.  No  special  antennas  or  spatial  processing  are  used  for  multipath  mitigation. 

6.  Doppler  effects  are  negligible. 

1.8  Approach 

This  thesis  presents  results  of  theoretical  analysis  and  simulation  of  multipath  mitigation 
using  the  complex  cepstrum  and  cepstral  domain  filtering.  Analyses  and  simulations  are  performed 
for  a  variety  of  multipath  scenarios  and  in  both  noiseless  and  additive  white  Gaussian  noise  en¬ 
vironments.  Computer  simulations  are  presented  in  Chapter  IV.  Simulations  are  written  using 
MATLAB  computational  software  from  The  Mathworks,  Inc.  of  Natick,  Massachusetts.  Simula¬ 
tions  are  performed  on  a  stand  alone  Pentium  based  PC  and  on  the  Sun  workstations  provided  by 
the  Air  Force  Institute  of  Technology  (AFIT). 
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II.  Cepstral  Analysis 


2.1  Introduction 

In  a  1963  paper  by  Bogert,  Healy,  and  Tukey,  entitled  “The  Quefrency  Alanysis  of  Time 
Series  for  Echoes:  Cepstrum,  Pseudoautocovariance,  Cross-Cepstrum,  and  Saphe  Cracking”  (1) 
it  was  observed  that  a  signal  composed  of  a  fundamental  and  echoes  of  the  fundamental  could 
be  decomposed  by  taking  the  logarithm  of  the  power  spectrum.  Bogert  et  ah  found  that  the 
logarithm  of  the  power  spectrum  contained  a  periodic  component  due  to  the  echoes  in  the  composite 
signal.  This  component  manifests  itself  as  periodic  peaks  when  the  Fourier  transform  is  taken. 
Because  their  technique  uses  the  spectrum  of  the  signal,  Bogert  et  al.  rearranged  the  letters  of 
spectrum  and  in  naming  their  new  technique  the  “cepstrum”.  Since  it’s  original  development, 
many  researchers  have  developed  ways  to  modify  the  cepstrum  for  application  to  specific  signal 
decomposition  problems.  To  avoid  confusion  with  these  new  techniques,  what  was  originally  called 
the  “cepstrum”  is  now  more  commonly  known  as  the  “power  cepstrum”.  Other  forms  of  the 
cepstrum  include  the  complex  cepstrum,  the  phase  cepstrum,  and  the  log  cepstrum.  This  thesis  will 
use  the  power  and  complex  cepstra.  As  a  final  note  of  introduction,  complex  cepstrum  techniques 
fall  into  a  class  of  nonlinear  filtering  techniques  developed  by  Oppenheim,  Schafer,  and  Stockham 
(7,  8)  termed  homomorphic  deconvolution. 

2.2  Definitions  of  the  Power  and  Complex  Cepstra 

2.2.1  Power  Cepstrum,  Definition.  As  previously  mentioned,  the  power  cepstrum  was  the 
first  cepstral  technique  developed.  This  technique  was  used  to  estimate  the  delay  and  amplitude  of 
echoes  present  in  a  composite  signal.  Since  computer  algorithms  are  typically  employed  to  evaluate 
the  power  cepstrum,  XpdnT],  it  is  usually  written  in  terms  of  the  Z-transform  of  a  sampled  signal 


11 


as  shown  in  Equation  8 


xv[nT]  =  {Z-^  {\og\X{z)\^)f  (8) 

where  X{z)  =  is  Z-transform  of  the  discrete  time  signal  x[nT],  n  is  the  sample 

number,  and  T  is  sampling  period.  The  logarithm  can  be  computed  relative  to  any  base;  however, 
the  natural  logarithm  is  most  commonly  used. 

Two  things  should  be  noted  concerning  the  power  cepstrum.  First,  if  |2;|  =  1  is  in  the  region 
of  convergence  of  X{z)^  then  the  counter  clockwise  contour  of  integration  is  typically  chosen  to 
be  the  unit  circle,  so  that  the  Z-transform  becomes  the  discrete  time  Fourier  transform  (DTFT). 
Unless  otherwise  noted,  the  DTFT  will  be  used  throughout  the  rest  of  this  thesis.  Second,  because 
the  logarithm  acts  only  on  the  magnitude  of  the  Z-transform,  the  phase  of  the  original  signal  is  lost. 
Therefore,  an  inverse  power  cepstrum  operation  cannot  be  applied  to  recover  the  original  signal. 
This  weakness  in  the  power  cepstrum  led  to  the  development  of  the  complex  cepstrum  which  will 
now  be  presented. 

2.2.2  Complex  Cepstrum  Definition.  The  complex  cepstrum  was  developed  out  of  the 
homomorphic  system  theory  of  Oppenheim,  Schafer,  and  Stockham  (7)  and  is  very  similar  to  the 
power  cepstrum.  The  primary  difference  is  that  the  complex  cepstrum  uses  the  complex  logarithm  of 
the  Z-transform,  maintaining  the  signal’s  phase  information.  The  general  definition  of  the  complex 
logarithm  and  complex  cepstrum,  x[nT]^  are  given  in  Equations  9  and  10,  respectively. 

logc  =  In  |c|  +  jlc  (9) 
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where  c  is  a  complex  number,  and  Zc  is  the  argument  of  the  complex  number,  c. 

f[nT]  =  <flog{X{z))z^-^dz  (10) 

^'^3  J  c 

where,  by  definition,  :r[0]  =  log(  .t[0]  ),  Finally,  as  a  matter  of  notation,  X{z)  will  be  used  to  denote 
logX(z). 

As  mentioned  earlier,  the  complex  cepstrum  is  more  versatile  due  to  the  fact  that  phase 
information  is  retained;  however,  carrying  this  phase  information  does  create  some  difficulty  in  the 
calculation  of  the  complex  cepstrum.  The  difficulty  arises  from  the  fact  that  the  complex  logarithm 
is  multi-valued  in  phase.  If  the  phase  of  the  logarithm  is  calculated  modulo  27r  (wherein  Zc  is 
the  principal  value  of  the  argument  of  c),  then  phase  discontinuities  result;  however,  this  is  not 
allowed  since  log[X(z)]  is  the  z-transform  of  x[nT]  and  thus  must  have  a  continuous  phase.  This 
difficulty  is  solved  using  any  of  several  readily  available  phase  unwrapping  techniques.  One  such 
technique,  presented  by  Childers,  Skinner,  and  Kemerait  (3)  adds  a  correction  term,  C{k)^  to  the 
phase  according  to  the  following  algorithm 

0  if  A;  =  0 

Cik  -  1)  -  27r  if  P{k)  -  P(^  -  1)  >  TT 
C{k)  (11) 

C{k  -  1)  H-  27r  if  P{k  -  1)  -  P{k)  >  tt 
C{k  —  1)  otherwise 

where  P{k)  is  the  phase  at  point  k. 

2,3  Cepstral  Processing  of  Multipath  Signals 

This  section  presents  a  typical  use  of  the  cepstrum,  namely  characterization  of  echoes  (such 
as  multipath)  in  a  composite  signal.  First,  the  theory  for  using  the  power  cepstrum  is  developed, 
followed  by  the  theory  of  the  complex  cepstrum.  Both  of  these  derivations  closely  follow  those  found 
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in  (2,  3).  Previous  work  has  shown  that  it  is  often  easier  to  detect  multipath  reflections  in  the  power 
cepstrum,  rather  than  the  complex  cepstrum  (3),  However,  the  complex  cepstrum  maintains  the 
needed  phase  information  contained  in  the  signal.  After  calculating  the  complex  cepstrum,  a  signal 
can  be  filtered  in  the  cepstrum  domain  to  remove  the  previously  detected  multipath  reflection, 
then  an  inverse  cepstral  operation  can  be  performed  to  recover  the  time  domain  direct  path  signal 
estimate.  This  process  will  be  described  in  more  detail  later. 

2.S.1  Power  Cepstrum.  To  see  how  the  power  cepstrum  can  be  used  to  decompose  a 
signal  comprised  of  a  direct  path  and  one  multipath  reflection,  let  x[nT]  be  the  signal  of  interest. 
This  signal  can  be  represented  as  follows 

x[nT]  =  f\nT]  *  g[nT]  (12) 

where  *  represents  the  convolution  operation,  f[nT]  represents  the  direct  path  signal,  and  g[nT]  is 
defined  as 

g\nT\  =  ^[nT]  +  aQ6[nT  —  tiqT]  (13) 

Referring  to  Equation  8  for  the  definition  of  the  power  cepstrum,  the  steps  in  calculating  the  power 
cepstrum  of  x[nT]  are  completed  in  the  following  manner 

X{z)  =  F{z)G{z)  (14) 

|X(^)|  =  \Fiz)\-\Giz)\ 

|X(^)p  =  |F(z)MG(z)P 

=  |F(^)P|(l+aoz-"«)P 
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Because  x[nT]  is  a  causal,  stable  signal,  the  DTFT  exists  and  is  found  by  substitution  of  z  = 


log|X(e^"^)P 


|F(e^“^)P  I  (1  +  ao  p 

log  |F(c^‘^^)p  +  log  I  (1  +  ao  P 

log  |F(e^"^)P  +  log  11  +  ao  cos(nowT)  -  jao  sin(now!r))P 

log  |JF(e^"^)P  +  log{[l  +  ao  cos(n.owT)]^  +  [ao  sin(noa;T)]^} 

log  |F(e^"^)P  +  log[l  +  2ao  cos(nowr)  +  a^  cos^(nowT)  +  a^  sin^(nowT)] 

log  |F(e^"^)P  +  log[l  +  2ao  cos(nowT)  +  al] 

log  |F(e^"^)P  +  log[(l  +  a^)(l  +  cos(nowT)] 

log  |F(c^"^)P  +  log(l  +  ao)  +  log[l  +  ^-2  cos(noa;T)] 

I  "h  Oq 

log  |F(e^-^)p  +  log(l  +  a")  +  log[l  +  +  e-^”o“^)]  (15) 


Next,  use  the  power  series  expansion,  log(l  x)  =  expand  the  third  term  in 

Equation  15.  Note  that  this  expansion  converges  only  if  |x|  <  1,  which  holds  for  this  application. 


To  prove  this,  let 


2ao 

1  +  ttn 


COS  (no  ci;  T) 


Then, 


0  ^  (1  —  ao  —  1  —  2  ao  +  aQ 


which  implies 


2  ao  1  4“  ®o 


or 


2ao 
1  -h  a^ 


<  1 
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Similarly,  (1  +  >  0?  leads  to  the  conclusion 


2  tto 
1  +  «o 


>  -1 


Therefore, 


-1  < 


2ao 


1  a? 


<  1 


(16) 


Because  |  cos  {no  ojT)\  <  1,  if  follows  that 


ki  =  I 


2  ao 
1  d-  ftn 


cos  (no  <^T)\  <  1 


(17) 


Therefore,  the  expansion  is  used  for  this  problem  and  is  given  by 


log 


^jnouT  ^  g 


-inowT^^  ^  ^  (-1) 


m+1 


m—1 


ao 


1  4-  ao 


^^jno(vT  _|_  g-~jnowT^ 


(18) 


Combining  Equations  15  and  18  gives 


log  |X(e^-^)p  =  log  |f’(e^"^)p  +  log(l  +  a^)  +  ^  +  e-^”''“^)]™  (19) 

m  i  “T  wq 

m—l 

=  log  |F(e^“^)p  +  log(l  +  al)  + 

1  +  tto 

-I  (i^y  (e'"""" +  5  (iTTl)’ 

=  log  |F(el“^)p  +  log(l  +  a^)  +  +  e-l"«“^) 

1  4-  ao 

_  1  / .gi2noa,T  ^-j2noC.T  3) 

2Vl  +  «o/ 

-(.1  I  _^2_  j  ^gjSnowT  _|_  ^-jZnQwT  _|_  ^^jnouT  ^^-jnouT'^  —  ,  . . 

3  \  1  +  ao  / 

Note  that  each  term  with  an  even  exponent  in  the  expansion  will  contain  a  constant,  with  decreasing 
magnitude  as  the  exponent  increases.  Also,  note  that  the  term  log(l  +  a^)  is  a  constant.  Lumping 
these  constants  together  into  another  constant  called  and  taking  the  inverse  Fourier  transform. 
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gives  the  power  cepstrum  of  this  signal. 


Xpa[nT]  =  +  + 

1  +  «0 


Uq 


2  V  1  + 


j2noUjT  _|_  ^~j2nowT 


+  2) 


f  ^  ^^jZnoi^T  ^-jSnoUjT  _|_  ^^jnoo^T  _|_  ^^-jnoujT^  —  •  •  •  } 

3  \  1  +  tto  / 

F-i{log  |F(e^“^)n  +  K6{t)  +  [5(<  +  noT)  +  6{t  -  noT)] 

i  +  ao 

+  2noT)  -  8{t  -  2noT)] 

z  1  -f  ag 

m  +  SnoT)  +  Sit  -  SnoT)]  -  . . . 


(20) 


Examination  of  Equation  20  shows  why  the  power  cepstrum  is  useful  for  determining  the  delay 
of  a  multipath  signal.  The  power  cepstrum  of  the  composite  signal  exhibits  periodic  peaks  which 
are  delayed  at  integer  multiples  of  the  multipath  delay;  thus,  the  power  cepstrum  can  be  used  to 
detect  a  multipath  component  in  a  composite  signal.  The  difficulty  occurs  when  one  attempts  to 
characterize  multiple  multipath  reflections.  This  difficulty  is  demonstrated  in  the  derivation  that 
follows,  which  parallels  a  derivation  in  (3). 

The  two  echo  case  is  derived  following  a  procedure  similar  to  that  for  a  single  echo.  In  this 
derivation,  arbitrarily  many  reflections  could  have  been  used.  However,  the  two  reflection  scenario 
is  useful  for  highlighting  the  key  points  without  adding  unnecessary  complexity. 

First,  refer  to  Equation  12  as  the  starting  point.  Then,  rewrite  Equation  13  to  account  for 
two  reflections,  which  yields  Equation  21 


g[nT]  =  S[nT]  -h  ao6[nT  -  noT]  +  aiS[nT  -  niT]  (21) 

where  and  n^T  are  the  amplitude  and  delay,  respectively,  of  the  reflection. 
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Next,  follow  the  same  steps  as  in  the  previous  derivation. 


X{z)  =  F{z)  •  G(z)  (22) 

\Xiz)\  =  |F(.)|.|G'(z)| 

|X(z)|2  =  \F{z)\^.\G{z)\-^ 

=  \Fiz)\^\{l  +  aoz-^«+aiz-^^)\^  (23) 


Evaluating  on  the  unit  circle,  z  =  e-’"*  and  taking  the  logarithm  gives 


|X(e^"^)|2 


=  I  (1  +  ao  +  ai  ^ 

=  log  \F{e^^^)\^  +  log  I  (1  +  ao  +  ai  ^ 


=  log|i^(6^’-^)P 


+  log  |1  -j-  ao  €08(710^^^)  +  ai  cos{niLjT)  —  j[ao  sin(noCt;T)  +  ai  sin(aia;T)]p 
=  \og\F{e^-'^)\^ 


+  log{[l  +  ao  cos(noa;T)  +  ai  cos(?^la;T)]^  +  [ao  sin(77oa;T)  +  ai  sin(7iiCL;T)]^} 
=  log  \F{e^^^)\^  +  log[l  H-  2ao  cos(noa?T)  +  2ai  cos(nia;T)  +  cos^(noa;T) 
+2aoai  cos(noa;T’)  cos(riia;T)  +  cos^{niUjT)  +  sin^(noCjT') 

+2aoai  sin(nowT')  sm(nia;T)  +  siii(niCjT)] 

=  log  \  F(e^^'^)\^  4-  log[l  4  ag  -t-  a^  -h  2ao  cos(aoCt;T)  4  2ai  cos(nia;T) 
+2aoai{cos[(ni  —  no)a;T]  4-  cos[(7ii  4-  no)6t;T] 

4-cos[(ni  -  no)a;T]  ~  cos[(ni  4  no)a;T]}] 


^  log,|i^(e^-^)P 


H-log{l  4  ag  4  0,1  4  2ao  cos(noa;T)  4-  2ai  cos{niUjT)  +  2aoai  cos[(ni  —  no)a;T]} 
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=  log  +  log  {  [1  +  ag  +  ai][l  +  ^  ^2 ^2  c«s(nowr) 

+1 + !§'+»?  1 + T+ »;  ~  * 

=  log|f'(e^"^)P  +  log(l  +  a^  +  a?) 

+ '”« I  ^ + i+»r+«t 


2tto^i 


1  +  +  ttj 


cos{(ni  “  ?^o)a;T} 


=  log  \F{e^‘^'^)f  +  log(l  +  +  a?)  +  log  {  1  +  j 


+  a^  +  a? 


j(eJnoa.T  ^g-ino«T) 

+ai(e^'”'"^  +  +  aoai(e^'^”'  ]  }  (24) 


Next,  the  final  term  of  Equation  24  can  be  expanded  in  a  similar  manner  as  for  the  single  refiection 
case.  Incorporating  the  expansion  into  Equation  24  gives 


F-i{log|Jf(e^“^)n  = 


F-qiog  f  {log(l  +  «o  +  «!)} 

=  F-Hlog|F(e^"^)n  +  J’-Mlog(l+«^+a2)} 


+F-^-, 


r[«o(e 


Jtiqu^T 


^  1  + 

2  1  +  fig  + 

+  i(- - J - 2  )'[ao(e''”°“^  + 

3  l-f  a5  4-  aj 

+aoai  -«o)a-T^]3  _ 

-  f -'{log  +  F-'  {  log(l  +  0.2  +  a?)} 

+f-'{ - ^ +  ai(e^”’‘"^  +  6“^”*“^) 

1+  ttg  + 
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2f«2/^i2nowr 


+  e 


-j2no<jjT 


+  2) 


+a2(ei2nia.T  ^  ^-j2n,^T  ^  3) 

+  e-^'2(ni-no)a.T  ^  2) 

_|_  g-i(’^i+^o)t^T  -no)a;T  _J_  g~i(ni-rio)a;T^ 

+a2ai(eJ"‘"^  +  -2«o)a-T  ^  -2no)u,T^ 


+aoa?(e^'"'>"^  +  g-inoo-T  ^  gi(2ni-no)a-T  ^  g-i(2ni -no VT^] 

+  i( - 1 - 2 

3  1  +  tio  “t" 


+  e-i3("i-«o)«T  ^  3gjnoa.T  ^  gg-inou-T^ 
+3o^ai(2e^'"’"^  +  +  e^'<2«o+niVT  +  g-i(2no+n,)«2’ 

+ei(2no-n,)<^T  ^  g-i(2no -nO<»T^  ^  3ao«? (2e^”<'"^  +  2e-^'"«“^ 
_|_gj(2ni+no)wT  _|_  g— j(2ni+no)wT  _|_  g;(2ni  — no)wT  _|_  g-j(2n] -no)wT^ 


+3ao®i(2e^^”'~”°^“^  +  +  e“J(3ni-no)a)r 

^ei(n,+no)o;r  ^  g-i(ni+no)wT^  ^  3aoai (2e^”""^  + 

_j_gi(3ni -2no)a;T  _j_  ^~j(3ni-2nQ)u;T  _|_  ^j(ni -2no)u;T  _|_  g-j(ni -2no)a;T ^ 
_|_^3^^(2gi(’^i— _|_  -|-  _|_  g-j(ni+no)wT 

_|_gi(ni -3no)a;T  g- j(ni -3no)wT^  _|_  +  ^j{2n-i -3no)u>T 

_|_^-j(2ni"no)a;r  _|_  gj(2ni -no)a;T  _|_  ^~j{2ni-no)(vT^  _|_  4^2^^ _j_  ^-j2nowT 
_^^j2niwT  _j_  g-j2nia;T  gj2(ni -no)a;T  _|_  j2(ni -no)a;T  _|_  2)]  _  .  _  j. 


(25) 
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Next,  the  inverse  Fourier  transforms  are  evaluated  to  yield  the  power  cepstrum.  After  collecting 
like  terms,  the  equation  above  becomes  the  following. 


(26) 


2  I  „2\  I  ^i)  ~l~  ^o)  +  (Xq  + 

1  ^2\3 


log(l  +  ttQ  -b  tti)  -j- 


(1  -\-al  +ai) 

6^0  +  18(1q  6uo  “b  (^^1  d~  13u|)  <io(3u^  -h  2IU2) 

6(1 +  ag  +a^)^ 

6^1  4-  I2a\  +  Gag  +  6ai  -  ao(3ai  -  6ai)  +  a0(6ai  +  ^a\) 

6(1  +  ttg  + 

— ci^(3ci^  —  12^1^  4"  3tt^  —  12fti)  H"  \2ciq(1\ 


m 

[S{t  —  uqT)  4-  S[t  4"  ^^0^)] 


6(1  4"  ttg  4"  ^1)^ 
balal  —  Soq  —  3aQ 


[<5(t  —  4“  ^{t  +  ^1^)] 


6(1  +  ^ +»;)>)  - 2”«^)  + «('  + 

3 

6a^ai  4-ar,('\2a 

4- 


4- 


+ 


3(l  +  a2  +  a2)3 
Gttgtti  4~  flo(12<xf  4"  17tti)  4“  “b  21a-j  4"  3tti) 
6(1  4-^0  +aiy 
ba^al  —  SaQai  —  3aoai 


[6[t  —  (jii  ~  no)T)  4'  S{t  4“  {ui  —  no)T’)] 


6(1  4-^0  4-af)3 


■  [^(i  —  2(rii  —  no)T)  -h  6(t  4"  2(ni  —  no)T)] 


alal 


[^(t  -  3(ni  -  no)r)  4-  S{t  4-  3(ni  -  no)T)] 


3(1  ^ao  4-  a^y 

+  — 2^1  2>t3  ~  ^0)^)  +  +  (^1  +  ^o)7")] 

d(1  4-  ag  4-  aiy 

_^2aoal  -  34ai -  3ag(af  +a0  ^ 

o(J-  r  CLq  4“  di) 

3ao(af  —  a?)  “  ^0^1 


4- 


4- 


4- 


4- 


4- 


6(1  4-  ^0  4" 


[5(i  —  (2ni  —  no)T)  4-  S(^t  +  {2ni  —  no)T)] 


a^ai 


(1  4-ag  +aiy 

al  ai 


[<5(i  —  (2no  4"  4-  ^{t  4-  (2no  4-  tii)T)] 


/I  .  2  ,  ~2W 

[1  -\-  ~r  (ii) 

[6{t  —  (2ni  4-  n{))T)  4-  +  (2ni  4-  nQ)T)] 


ao«i 


(1+al+aiy 

apal 

(1  -\-a\  +af)® 


[^(i  —  (3tii  —  77-1)7’)  +  6{f,  +  (3771  —  770)7")] 
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+  -  2^0)^)  +  +  (3^1  -  2»^o)r)] 

v-*-  >  ®o  ^1) 

+  3(l+1fiaf)3 

3  2 

+  ls(t  -  i2m  -  3no)r)  +  S{t  +  (2m  -  3no}T)]  +  ... 

»3^I  “T  ^0  '  ^1) 

From  Equation  26  it  is  seen  that  the  power  cepstrum  will  contain  peaks  at  the  delays  of  each 
reflection,  just  as  it  did  in  the  single  reflection  case.  However,  it  is  also  seen  that  the  reflections 
will  interact  with  each  other  causing  peaks  at  sums  and  differences  of  the  delays.  It  is  these  peaks 
that  cause  the  difficulty  in  determining  the  true  multipath  delays.  If  a  sufficiently  small  number  of 
multipath  reflections  exist,  the  true  delays  can  be  distinguished  from  these  interference  delays  using 
an  exhaustive  search.  This  method  would  involve  comparing  the  delay  associated  with  every  peak 
to  the  sums  and  differences  of  the  other  delays,  and  determining  which  delays  are  true  multipath 
delays.  Using  a  computer  algorithm,  this  procedure  could  probably  be  implemented  for  cases 
where  only  a  few  reflections  exist.  However,  this  method  would  become  very  time  consuming  if 
many  reflections  exist.  Therefore,  rather  than  dealing  with  this  difficulty,  this  thesis  will  focus  on 
using  the  complex  cepstrum  to  remove  all  reflections  without  necessarily  estimating  the  parameters 
associated  with  the  reflections. 

2.3.2  Power  Cepstrum  Example.  To  see  how  the  power  cepstrum  can  be  used  to  determine 
the  delay  in  a  single  multipath  reflection  case,  consider  the  following  simple  example.  Let  x[nT]  be 
the  direct  path  signal 

x[nT]  =  (27) 

where  T  is  the  sampling  period,  n  is  the  sample  number,  and  the  signal  is  given  an  initial  amplitude 
of  1  so  that  everything  is  normalized  to  the  direct  path. 

Next,  let  a  single  multipath  reflection  occur  three  time  units  later  than  the  direct  path,  with 
an  amplitude  of  one  half  with  respect  to  the  direct  path.  Then,  from  Equation  6  the  composite 
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signal  is  given  by 


c[nT]  = 


(28) 


where  c[nT]  is  the  composite  signal.  A  plot  of  these  signals  is  shown  in  Figure  3. 


Composite  Signal 


Figure  3  Signals  for  Power  Cepstrum  Example 


Following  the  procedures  previously  developed,  the  power  cepstrum  of  the  composite  signal 
was  calculated  and  is  shown  in  Figure  4.  It  is  seen  in  Figure  4  that  the  power  cepstrum  contains 
peaks  at  ±3  time  units,  as  predicted  by  Equation  20. 


2.S.3  Complex  Cepstrum.  The  development  of  the  complex  cepstrum  in  a  multipath 
environment  is  very  similar  to  that  of  the  power  cepstrum.  The  difference  lies  in  the  fact  that  the 
complex  logarithm  is  used  to  retain  phase  information.  As  mentioned  earlier,  this  allows  one  to 
filter  the  signal  in  the  cepstral  domain,  then  perform  an  inverse  operation  to  return  to  the  time 
domain.  This  j)rocess  is  shown  in  Figure  5.  The  single  reflection  complex  cepstrum  derivation 
follows.  As  before,  this  derivation  parallels  those  in  (2,  3). 
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Rea!  Cepstrum  of  Composite  Signal 


Time 


Figure  4  Power  Cepstrum  for  Example  Signal 


Figure  5  Complex  Cepstrum  Filtering  Process 


For  the  complex  cepstrum  derivation,  begin  with  a  signal  which  is  the  convolution  of  two 
other  signals. 


x[nT]  =  f[nT]  *  g[nT] 


(29) 


where  *  represents  the  convolution  operation,  f[nT]  represents  the  direct  path  signal,  and,  for  a 
single  reflection,  g[nT]  is  defined  as 


g[nT]  =  6[nT]  +  aoS[nT  —  uqT] 


(30) 
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Referring  to  Equation  10,  the  steps  in  computing  the  complex  cepstrum  are  completed  as  follows. 


Xiz)  =  Fiz)G{z) 

=  F{z){l  +  aoz-^^°) 


(31) 


Evaluating  the  Z-transform  on  the  unit  circle  and  taking  the  logarithm  yields 

(32) 

(33) 

Again,  a  power  series  expansion  is  used  for  the  final  term  above,  provided  that  |ao  = 

|ao|  <  1,  which  requires  that  reflections  be  attenuated  relative  to  the  direct  path  signal.  This 
generally  is  true  due  to  the  conservation  of  energy,  but  can  be  violated  when  many  reflections 
combine  in  phase  to  form  a  composite  reflection  of  greater  magnitude. 

.0.1(1+ (34) 

k=l 

Using  the  expansion  yields 

log[X(e^-^)]  =  log[F(e^'"^)]  +  f;(-l)^+^  (35) 

k=l 

Finally,  the  inverse  Fourier  transform  is  taken  to  give  the  complex  cepstrum. 

i^-i{log[X(e>^)]}  =  i^-i{log[F(e^-^)]} 

-  . . .  |  (36) 

2  3  J 

2  3 

=  f -i{log[F(e^“^)]}  +  aoS{t  -  noT)  -  "^6(1  -  2noT)  +  ^6{t  -  3noT)  -  . . . 


X(e3<^^)  =  F(e3^)(l  +  aoe-3“”o^) 

log[X(e3“^)]  =  log[F(e3-’’)]  +  log[(l  +  aoe-^‘^^'^)] 
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Just  as  in  the  power  cepstrum  case,  it  is  seen  that  the  complex  cepstrum  will  contain  peaks  at 
integer  multiples  of  the  multipath  delay.  However,  unlike  with  the  power  cepstrum,  these  peaks 
can  be  removed  from  the  complex  cepstrum  by  filtering.  Following  this  filtering  operation  with  an 
inverse  complex  cepstrum  will  then  yield  a  multipath  free  signal  in  the  time  domain. 

The  same  procedure  as  above  can  be  followed  for  the  multiple  echo  scenario  (3).  Again,  as  in 
the  power  cepstrum  case,  a  two  echo  case  will  be  presented  for  illustration  purposes.  However,  the 
analysis  could  be  extended  to  the  case  of  arbitrarily  many  echoes. 

Again,  start  with  a  signal  which  is  the  convolution  of  two  signals,  as  given  in  Equation  29. 
However,  in  this  case  define  g[nT\  as 

g[nT]  =  d[nT]  +  ao^[n-T  —  uqT]  +  ai8[nT  —  n^T]  (37) 

where  ai  and  UiT  are  the  amplitude  and  delay,  respectively,  of  the  multipath  signal. 

Then,  the  derivation  follows  as  before. 


X(z)  =  Fiz)Giz) 

=  F{z){l  +  aoz-'^<‘  +aiz-'^') 


(38) 


Evaluating  the  Z- transforms  on  the  unit  circle  and  taking  the  logarithm  yields 


log[X(e>^)]  =  log[F(e^'"^)]  +  log[(l  +  aoe-'""«^  +  aie-^"”^'^)] 


(39) 


Using  the  power  series  expansion  for  the  last  term,  when  |  ao  e  +  ai  e  |  <  i.  gives  the 

following 


log[X(e^“^)]  =  log[i^(e^“^)]  + 

k=l 


k 


(40) 
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Finally,  taking  the  inverse  Fourier  transform  yields  the  complex  cepstrum. 


F-Mlog[X(e^"^)]} 


F-i{log[F(e^“^)]} 

+F-^  {  ^ ^ ^ 

+  3  •••} 

f -illoglFCe^'""^)]}  +  F-^  {  aoe-^“”o^  +  -  iage 

+  ia?e 

2  2  3  u 

+  aoa?e-^'"("'>+2ni)T  ^ 


(41) 

i2 


—jcjnoT 

—jSivn^T 


F  ^{log[F(e^"^)]}  +  aoS{t  -  rioT)  +  ai(5(i  -  riiT) 

—  ^al6{t  —  2noT)  -  ^a\6{t  -  2n\T)  —  a\ai6{t  —  (no  +  n{)T) 

-\-^alS(t  ~  SnoT)  +  —  3niT)  +  alai6{t  ~  (2no  +  ni)T) 

o  o 

+aoOi^(i  “  (^^0  "1“  271i)T)  —  . . . 


As  in  the  single  reflection  case,  it  is  seen  from  Equation  41  that  the  complex  cepstrum  will  contain 
peaks  with  amplitudes  and  delays  proportional  to  those  of  the  multipath  signals.  Also,  like  the 
power  cepstrum,  it  is  seen  that  the  multiple  reflections  will  interact  causing  difficulty  in  using  the 
complex  cepstrum  to  characterize  the  multipath  signals.  However,  unlike  in  the  power  cepstrum 
case,  the  complex  cepstrum  can  be  filtered  in  the  cepstral  domain,  then  inverted  back  to  the 
time  domain.  Using  this  technique,  the  multipath  signals  are  removed  and  the  direct  path  signal 
recovered  without  needing  to  fully  characterize  the  multipath  signals. 

2.4  Filtering  in  the  Cepstral  Domain 

Three  methods  for  filtering  a  cepstral  domain  signal  are  long  pass,  short  pass,  and  comb 
filtering.  Each  of  these  methods  has  advantages  and  disadvantages  depending  upon  the  goals  of 
the  designer.  The  long  pass  filter  is  the  cepstral  equivalent  to  a  frequency  domain  high  pass  filter. 
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The  long  pass  filter  is  designed  to  set  all  points  in  the  cepstrum  prior  to  the  first  echo  peak  equal  to 
zero.  Thus  using  a  long  pass  filter  preserves  the  multipath  signals  while  filtering  out  the  direct  path 
signal.  Similarly,  the  short  pass  filter  is  the  cepstral  equivalent  of  a  low  pass  filter.  The  short  pass 
filter  replaces  everything  from  the  first  echo  peak  forward  (in  time)  with  a  zero.  Thus,  the  short 
pass  filter  removes  all  multipath  signals  while  passing  the  direct  path  vsignal.  The  disadvantage 
of  both  the  long  pass  and  short  pass  filters  is  that,  by  zeroing  out  large  portions  of  the  complex 
cepstrum,  these  filters  may  cause  significant  degradation  to  the  signal (s)  of  interest.  An  alternative 
to  these  methods  is  the  comb  filter.  This  filter  is  the  cepstral  equivalent  to  a  notch  filter.  The 
comb  filter  replaces  multipath  peaks  in  the  complex  cepstrum  with  the  average  of  the  two  (or 
more)  points  immediately  adjacent  to  the  peaks.  In  this  manner,  the  comb  filter  removes  the  delta 
functions  in  the  complex  cepstrum  which  are  caused  by  the  multipath  interference.  Thus,  when 
the  filtered  signal  is  converted  back  to  the  time  domain,  the  direct  path  signal  is  recovered.  The 
difficulty  with  this  technique  lies  in  detecting  the  peaks  in  the  complex  cepstrum.  As  previously 
mentioned,  researchers  (3)  have  found  that  it  is  often  easier  to  detect  multipath  peaks  by  using  the 
power  cepstrum  rather  than  the  complex  cepstrum.  A  parallel  process  is  employed  to  detect  and 
mitigate  multipath  effects;  the  power  cepstrum  is  used  to  detect  multipath  delays  and  the  complex 
cepstrum  is  used  to  filter  out  the  multipath  reflection.  A  block  diagram  of  this  process  is  shown  in 
Figure  6,  A  simple  example  using  the  complex  cepstrum  to  filter  out  multipath  follows. 

2.4-1  Multipath  Mitigation  Example.  As  an  example  of  the  theory  that  has  been  presented 
consider  again  the  decaying  exponential  signal,  x[nT] 

x[nT]  =  (42) 

where  T  is  the  sampling  period,  n  is  the  sample  number,  and  the  signal  was  given  an  initial 
amplitude  of  1  so  that  everything  is  normalized  to  this  signal. 
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Figure  6  Power  and  Complex  Cepstrum  Multipath  Detection  and  Mitigation  Procedure 


If  x[nT]  is  the  direct  path  signal,  and  a  multipath  reflection  occurs  three  time  units  later  with 
an  amplitude  of  one  half  with  respect  to  the  direct  path,  then  the  composite  signal  is  represented 
as  follows. 


c[nT]  =  x[nT]  +  0.5a;[nT-3]  (43) 

c\nT\  = 

where  c[nT]  is  the  composite  signal.  A  plot  of  these  signals  is  shown  in  Figure  7. 

Recall  that  the  procedure  for  removing  the  multipath  signal  is  shown  in  Figure  5.  The  first  step  is 
to  calculate  the  complex  cepstrum  of  the  composite  signal  as  previously  presented.  A  plot  of  the 
complex  cepstrum  for  this  example  is  shown  in  Figure  8.  Note  that  the  complex  cepstrum  contains 
peaks  at  3,  6,  and  9  time  units  as  predicted  by  the  theory  previously  presented.  The  next  step 
in  removing  the  multipath  signal  is  to  comb  filter  the  complex  cepstrum  at  time  units  3,  6,  and 
9.  This  is  accomplished  by  replacing  complex  cepstrum  values  at  those  times  with  the  averages  of 
the  two  points  on  either  side  of  the  peaks.  The  complex  cepstrum  of  the  multipath  signal  and  the 
comb  filtered  complex  cepstrum  are  both  shown  in  Figure  8.  Finally,  the  comb  filtered  complex 
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Composite  Signal 


Figure  7  Signals  for  Complex  Cepstrum  Example 


cepstrum  is  transformed  back  into  the  time  domain  using  an  inverse  complex  cepstrum  operation 
to  recover  the  direct  path  signal.  The  recovered  direct  path  signal  for  this  example  is  shown  in 
Figure  9  along  with  a  plot  of  the  original  direct  path  signal.  The  mean  squared  error  between  the 
recovered  direct  path  signal  and  the  actual  direct  path  signal  is  1.16  x  10 demonstrating  that 
this  method  is  effective  for  multipath  mitigation  for  this  signal. 
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IIL  Complex  Cepstrum  Adaptive  Filter 


3. 1  Overview 

This  chapter  describes  the  Complex  Cepstrum  Adaptive  Filter  (CCAF)  and  the  non-coherent 
delay  lock  loop,  also  commonly  called  an  early-late  gate.  Additionally,  the  chapter  describes  the 
input  signal  model.  The  first  section  covers  the  signal  model.  The  subsequent  sections  discuss  the 
CCAF  and  it’s  components.  The  final  section  describes  the  NCDLL. 

3.2  Received  Signal  Model 

Prior  to  processing  the  received  signal  for  multipath  mitigation,  it  is  assumed  that  the  signal 
has  been  mixed  down  to  some  convenient  intermediate  frequency  (IF).  If  the  received  multipath 
signal  is  assumed  to  have  one  refiection,  the  received  signal  can  be  described  by  the  following 
equation: 

s{t)  =  -  Td)  cos  (ojijrt  -f-  ^o)  +  \/2Pai  c{t  —  Tm)  cos  {iViFt  +  ^i)  +  n(t)  (44) 

where  ao,  ai  are  the  direct  path  and  multipath  amplitudes,  respectively,  Oq  =  —cvoTd^  0i  = 
are  the  direct  path  and  multipath  phases,  respectively,  and  Tm  are  the  direct  path  and  multipath 
propagation  delays,  respectively,  cjq  is  the  carrier  angular  frequency,  and  n{t)  is  additive  white 
Gaussian  noise  (AWGN).  Both  phases  are  measured  relative  to  the  phase  of  the  transmitter  at  the 
time  of  transmission,  which  is  assumed  to  be  zero  without  loss  of  generality. 

3.3  Complex  Cepstrum  Adaptive  Filter 

3.3.1  Overview.  A  block  diagram  of  the  CCAF  is  shown  in  Figure  10.  This  unit  uses 
the  complex  cepstrum  in  parallel  with  the  power  cepstrum  as  described  in  Chapter  II.  The  system 
calculates  the  power  cepstrum  for  detection  of  the  peaks  associated  with  the  multipath  signal.  The 
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peak  delays  are  then  provided  to  the  adaptive  weighted  comb  filter  which  removes  the  corresponding 
peaks  from  the  complex  cepstrum.  An  inverse  complex  cepstriim  is  calculated,  yielding  a  replica  of 
the  direct  path  signal.  The  fidelity  of  the  direct  path  replica  depends  on  several  factors  including 
signal-to-noise  ratio  and  aliasing  effects.  These  factors  will  be  discussed  in  the  sections  that  follow. 
It  is  proposed  that  the  CCAF  be  inserted  into  a  typical  GPS  receiver  code  tracking  loop  as  shown 
by  the  dotted  box  in  Figure  10.  By  placing  the  CCAF  in  this  position,  multipath  interference  can 
be  removed  from  the  composite  signal  prior  to  code  tracking,  thus  eliminating  code  tracking  error 
caused  by  the  multipath. 


Figure  10  A  Code  Tracking  Loop  Employing  the  Complex  Cepstrum  Adaptive  Filter 

3.S.2  Zero  Padding.  The  zero  padding  block  in  Figure  10  inputs  a  block  of  data  (taken 
to  be  one  period  of  code  for  this  work),  and  appends  zeros  to  the  end  of  the  data  record.  This 
has  two  effects  on  the  data  processing.  First,  if  the  correct  number  of  zeros  are  appended,  the 
record  length  can  be  increased  to  be  a  power  of  two.  This  allows  computation  of  the  discrete 
Fourier  transform  (DFT)  and  the  inverse  DFT  using  faster  computation  algorithms,  thus  speeding 
up  processing  time.  Second,  as  discussed  in  (3),  appending  zeros  increases  the  frequency  domain 
resolution.  This  increase  reduces  or  eliminates  cepstral  (pseudo-time)  domain  aliasing,  and  reduces 
phase  unwrapping  errors. 
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3.3.3  Cepstra  Calculations.  Following  zero  padding,  the  power  and  complex  cepstra  of 
the  data  record  are  calculated.  As  mentioned  earlier,  both  cepstra  are  calculated  because  it  is 
often  easier  to  detect  multipath  peaks  using  the  power  cepstrum  rather  than  the  complex  cepstrum 
(3);  however,  the  complex  cepstrum  must  be  used  for  time  domain  signal  reconstruction  following 
filtering.  These  calculations  are  performed  using  the  definitions  given  in  Equations  8  and  10.  The 
calculations  are  covered  extensively  in  Chapter  II,  and  will  not  be  discussed  again  here. 

3.3.4  Peak  Detection.  In  order  to  remove  the  reflection  from  the  complex  cepstrum,  the 
first  peak  associated  with  the  reflected  signal  must  be  detected.  A  magnitude  plot  of  a  typical 
power  cepstrum  for  a  GPS  multipath  signal  is  shown  in  Figure  11.  The  first  delta  function  in  the 
power  cepstrum  corresponds  to  the  delay  of  the  reflected  signal,  assuming  a  single  reflection,  and 
provides  all  the  information  necessary  to  filter  out  the  reflection.  The  peak  detector,  of  Figure  10, 
performs  the  function  of  determining  the  delay  associated  with  the  first  peak.  In  Chapter  II,  it  is 
shown  that  the  power  and  complex  cepstra  contain  delta  functions  at  the  multipath  delay  and  all 
integer  multiples  of  the  delay.  Therefore,  the  peak  detector  must  be  able  to  detect  the  first  delta 
function  for  the  single  reflection  case.  To  perform  this  task,  the  peak  detector  compares  the  power 
cepstrum  value  at  each  point  to  the  two  adjacent  points.  A  point  which  has  a  magnitude  exceeding 
those  of  the  two  adjacent  points  by  more  than  some  user  set  threshold  is  declared  a  multipath 
peak.  In  Equation  20,  it  is  seen  that  the  first  delta  function  in  the  power  cepstrum  will  have 
magnitude  |ao/  (l  +  ^o)  1^1/2  (see  Equation  16)  and  each  subsequent  delta  will  have  magnitude 
[ao/  (1  +  ao)]*"  where  k  is  an  integer  corresponding  to  the  harmonic  number  of  the  delta  function 
(ie.  k  —  1,  2,  3  . . .).  Therefore,  the  first  peak  will  have  a  greater  magnitude  than  each  subsequent 
peak.  The  peak  detector  searches  for  the  first  peak  present  and  declares  the  corresponding  delay 
to  be  that  of  the  reflection.  As  shown  in  Figure  11,  the  power  cepstrum  will  contain  peaks  at  the 
chip  interval  as  well  as  the  multipath  delay.  However,  as  discussed  in  the  next  section,  the  adaptive 
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weighted  comb  filter  is  designed  not  to  filter  these  peaks.  Therefore,  the  detector  is  designed  not 
to  detect  a  delay  of  exactly  1,0  chip. 


Magnitude  Plot  of  the  Power  Cepstrum 


Figure  11  Typical  Power  Cepstrum  Magnitude  for  a  GPS  Multipath  Composite  Signal 

S.S.5  Adaptive  Weighted  Comb  Filter.  The  adaptive  weighted  comb  filter  performs  the 
function  of  removing  the  reflected  signal  from  the  complex  cepstrum.  Following  detection  by  the 
peak  detector,  the  delay  of  the  reflection  is  provided  to  the  adaptive  weighted  comb  filter.  This 
filter  then  replaces  the  complex  cepstrum  data  point  at  the  detected  delay  with  the  average  of  the 
previous  point  and  the  next  point.  Additionally,  the  comb  filter  replaces  the  data  points  at  each 
integer  multiple  of  the  multipath  delay  with  the  average  of  the  two  adjacent  points.  In  this  manner 
the  delta  functions  in  the  complex  cepstrum,  due  to  the  reflected  signal,  are  x'eplaced  by  linearly 
interpolated  estimates  of  the  complex  cepstrum  values  for  the  direct  path  signal.  A  portion  of  a 
typical  filtered  and  unfiltered  complex  cepstrum  is  shown  in  Figure  12.  Observe  in  Figure  12  that 
the  complex  cepstrum  contains  delta  functions  at  the  chip  interval,  as  well  as  those  due  to  the 
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multipath  signal.  Simulations  show  that  these  peaks  are  important  for  proper  reconstruction  of 
the  direct  path  estimate.  Therefore,  the  adaptive  weighted  comb  filter  is  designed  not  to  filter  the 
complex  cepstrum  at  the  chip  interval. 


Filtered  and  TJiifiltered  Complex  Cepstrum 


Figure  12  Typical  Filtered  and  Unfiltered  Complex  Cepstrum  for  a  GPS  Multipath  Composite 
Signal 

To  better  understand  the  adaptive  weighted  comb  filter,  consider  a  complex  cepstrum  data 
sequence,  f  (nT),  where  n  is  the  sample  number  and  T  is  the  sampling  interval.  Assume  that  x{nT) 
is  N  samples  long,  and  contains  an  initial  multipath  peak  in  the  Ck  element.  Then,  x{nT)  could  be 
written  as  follows: 


x{nT)  =  [  Cl  C2  ...  Cfc_i  Ck  Cfc+l  .  .  .  C2k-1  C2k  C2k-\-l  •  •  •  CJV  ] 


where  Ci  is  the  data  point  in  the  complex  cepstrum. 
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Next,  a  filter  vector  W,  also  of  length  N,  is  defined  as  follows: 

Cfc-l  H-  Cfc+I  ^  ^  C2k~l  +  C2k+1  ^  ^ 

2ca;  2c2fc 

where  W  consists  of  ones  for  all  elements  except  those  corresponding  to  the  multipath  peaks  in 
x{nT)  (i.  e.  the  elements  k,  2k,  3k,  etc.).  It  is  easily  shown  that  if  the  vectors  x{nT)  and  W 
are  multiplied  in  a  point-wise  fashion,  the  result  is  a  filtered  complex  cepstrum  vector  with  the 
multipath  reflections  removed. 

S.3.6  Inverse  Complex  Cepstrum.  After  removing  the  multipath  reflection,  a  time  domain 
estimate  of  the  direct  path  signal  is  recovered  through  an  inverse  complex  cepstrum  operation.  The 
steps  of  this  calculation  are  shown  in  Figure  5  of  Chapter  II.  The  process  is  simply  the  inverse  of  the 
complex  cepstrum  process.  First,  a  forward  Fourier  transform  is  calculated,  followed  by  a  complex 
exponentiation  (with  phase  wrapping),  and  finally  an  inverse  Fourier  transform  is  calculated  yielding 
a  time  domain  signal.  The  phase  wrapping  procedure  is  the  inverse  of  the  chosen  phase  unwrapping 
procedure  used  in  the  complex  cepstrum  calculation,  as  discussed  in  Section  2.2.2. 

S.4  The  Non- Coherent  Delay  Lock  Loop 

Tracking  of  the  GPS  spreading  code  is  a  primary  function  of  the  GPS  receiver.  This  tracking 
provides  a  method  of  determining  the  line-of-sight  distance  from  the  receiver  to  the  satellite,  and 
also  acts  to  despread  the  received  signal.  Code  tracking  is  typically  accomplished  using  a  delay 
lock  loop  (DLL),  frequently  called  an  early-late  gate.  The  delay  lock  loop  may  operate  in  either 
a  coherent  or  non-coherent  fashion.  Since  coherent  DLLs  require  an  estimate  of  carrier  phase, 
considerable  interaction  between  the  carrier  and  code  tracking  loops  is  necessary  for  a  coherent 
DLL  to  function  properly.  Additionally,  cycle  slii)s  in  the  carrier  tracking  loop  can  cause  a  loss  of 
code  tracking  lock  when  a  coherent  DLL  is  used.  For  these  reasons,  the  coherent  DLL  is  considered 
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somewhat  fragile  and  typically  not  used  for  GPS  applications  (9).  Therefore,  the  non-coherent 
delay  lock  loop  will  be  presented  here.  For  analysis  purposes,  a  block  diagram  of  a  simplified 


S(iiiaro  Law  Detector 


Figure  13  The  Non-Coherent  Delay  Lock  Loop 


NGDLL  is  shown  in  Figure  13.  In  this  DLL,  the  received  signal  is  cross  correlated  with  early  and 
late  versions  of  the  locally  generated  spreading  code  replica.  The  results  of  these  correlations  are 
then  bandpass  filtered,  squared,  lowpass  filtered,  and  subtracted  to  form  the  discriminator  output 
for  the  loop.  To  close  the  tracking  loop,  the  discriminator  output  is  filtered  via  a  loop  filter,  then 
input  to  a  voltage  controlled  clock  which  triggers  the  PN  sequence  generator  producing  the  on  time 
code  replica. 

3. 4-0.1  Signal  Component  NCDLL  Analysis.  The  following  analysis  patterns  the 
development  of  (10).  Consider  a  DS/SS  signal 

r(t)  =  Vweit  -  Td)  cos[wo<  +  0d{t  -  Td)  +  ^]  +  n(t)  (45) 
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where  r{t)  is  the  received  signal,  P  is  the  received  power,  is  the  signal  propagation  delay, 
c{t  -  Td)  is  the  delayed  spreading  code,  ojq  is  the  carrier  radian  frequency,  Od{t  —  is  the  delayed 
data  signal,  <f>  is  an  arbitrary  carrier  phase  angle,  and  n{t)  is  bandlimited  zero  mean  white  Gaussian 
noise,  introduced  by  the  channel. 

Assuming  a  two  sided  noise  power  spectral  density  of  7Vo/2  W /Hz,  the  noise  can  be  represented 
in  terms  of  its  in-phase,  nj{t)^  and  quadrature,  nQ{t)^  components. 

n{t)  =  ^/2nr{t)  cos{(jJot)  -  V2nQ{t)  sm{(jJot)  (46) 

After  power  division,  the  received  signal  in  each  correlator  arm  becomes 

P(t)  =  VPc{t-rd)  cos[a;ot  +  Odit~Td)  (f>]  ni{t)  cos{u;ot)  -  7iQ{t)  sin(a;ot)  (47) 

This  signal  is  subsequently  correlated  with  early  and  late  replicas  of  the  spreading  code.  Assuming 
that  the  NCDLL  operates  at  a  fixed  intermediate  frequency,  the  early  and  late  locally  generated 
signals  are  given  by  Equations  48 

asit)  =  ^/Klc^  -  fd  +  Y^c)  cos[(ct>o  -  (^iF)t  +  <p']  (48) 

aiit)  =  (i  -  Td  -  ^^<=)  <^os[(a;o  -  ujif)*  +  4>'] 

where  and  ai  are  the  early  and  late  correlation  signals,  respectively,  Ki  is  an  RF-to~IF  conversion 
loss  constant,  fd  is  the  loop’s  estimate  of  the  propagation  delay,  A  is  the  correlator  spacing  in  chips, 
and  lojf  is  an  intermediate  frequency. 
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The  early  correlator  output  is 


yeit)  =  r'{t)  ■  asit)  (49) 

-  -  Td)  cos  [wo^  +  Od{t  -  Td)  +  <!>]+  ni{t)  cos(woi)  -  sin(wot)j  • 

c(^  -  Td  y cos  [(wo  -  ^if)^  +  <t>'\ 

Using  appropriate  trigonometric  identities,  the  above  equation  may  be  simplified  by  noting  that 
yE  {t)  is  filtered  by  a  bandpass  filter  centered  at  the  intermediate  frequency  to  obtain  ze  (t).  There¬ 
fore,  only  the  difference  terms  of  the  simplification  need  be  considered.  Thus,  the  early  correlator 
output,  ze  (t),  becomes 


Zsit) 


y/PKi  c{t  -  Td)  c{t  -Td+  ^Tc)  cos  (wjj-i  +  0d{t  -  Td)  +  4>  -  </)' )  (50) 

+  ^/kI c{t  -  Td  +  y  Tc)  \ni{t)  cos  {wipt  -  4>')  -  nq{t)  sin  {ujipt  -  0')] 


Using  a  similar  process,  the  late  correlator  signal,  zi{t)^  may  be  expressed  as 


ZL{t) 


\/PKi  c{t  -  Td)  c(t  -Td  -  y  ^c)  cos  (wipt  +  0d{t  -Td)  +  (f>  -  ^')  (51) 

+  c(t  -Td  -  yTc)  [ni(t)  cos  {ivipt  -  0')  -  nq(t)  sin  {u/pt  -  4>')] 


From  the  previous  two  equations,  it  seen  that  the  signals  ze  (i)  and  zl  {t)  are  composed  of 
the  desired  signal  and  the  noise  terms,  including  code  self  noise^.  If  the  spread  spectrum  processing 
gain  is  sufficiently  high  (typically  greater  than  10),  the  self  noise  term  can  be  neglected.  According 
to  (10),  “The  amount  of  performance  improvement  that  is  achieved  through  the  use  of  spread 
spectrum  is  defined  as  the  processing  gain  of  the  spread  spectrum  system.’'  The  processing  gain  is 

^Self  noise  is  defined  as  sn  —  E{c{t  +  T)c'{t  -|-  f)]  —  c(t  +  r)c'{t  +  f)  where  E  is  the  expected  value,  c{t  +  t) 
is  the  spreading  code,  c'(t  -|-  f)  is  the  derivative  of  the  spreading  code  with  respect  to  time,  r  is  the  propagation 
delay,  and  f  is  the  NCDLL  estimate  of  the  propagation  delay.  Self  noise  is  a  broadband  noise-like  waveform  (9). 
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typically  defined  to  be  the  spread  spectrum  bandwidth  divided  by  the  data  bandwidth.  Using  this 
definition,  the  processing  gain  for  GPS  C/A  code  can  be  calculated  as  follows 

—  50 

“  1 

1.023  X  10^ 

=  20460 

where  T  is  the  data  bit  time,  and  Tc  is  the  chip  width,  both  in  seconds. 

From  Equation  52,  it  is  seen  that  the  processing  gain  for  GPS  C/A  code  is  indeed  high,  and 
therefore  the  code  self  noise  component  can  be  neglected.  Ignoring  the  self  noise  component  and 
the  AWGN,  the  noise-free  components  of  ze  and  zl  are 

=  VKPc{t  -  rd)c{t  -  Td  +  Tc)  cos  [wjirt  +  0d(t  -  Td)  -f  -  <j)']  (53) 

ZsL  (t)  =  \/ATc(t  -  Td)c{t  -  Td  ~  Tc)  COs[cJiFt  +  0d(t  -  T^)  -  0'] 

where  the  “s”  notation  is  used  to  distinguish  from  the  more  complete  signals  of  Equations  50  and 
51. 

The  dc  component  of  the  spreading  code  multiplication,  by  definition,  is  the  code  autocorre¬ 
lation,  Rc{t),  evaluated  at  r  =  -f  (A/2)  Tc  in  the  case  of  and  r  =  Td  —  Td  —  (A/2)  Tc 

in  the  case  of  zl]  Td  is  the  propagation  delay,  fd  is  an  estimate  of  the  propagation  delay,  and  A  is 
the  correlator  spacing  in  chips.  Defining  the  NCDLL  tracking  error  as  ^  =  (r^  —  fd)  /Tc^ 

ZsEit)  =  VKPRc  (5  +  I)  Tc  cos  [uJiFt  +  Odit  -Td)  +  <p-  <p']  (54) 

ZsL{t)  =  VKPRc  (<^  -  y)  Tc  COs[(JiFt  +  0d{t  -  Td)  +  4>  -  (p'] 
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Assuming  the  bandpass  filters  are  designed  to  pass  IF  signals  with  no  distortion  of  the  messages, 
the  input  to  the  square-law  envelope  detectors  is  given  by  Equation  54  above.  Also,  assuming 
that  the  low  pass  filters  completely  block  signal  components  at  signal  component  of  the 

NCDLL  discriminator  can  be  shown  to  be  (10) 

e{t,6)  =  [xlit)  -  xlit)]  (55) 

=  5 I)  +y)^} 

Defining  the  NCDLL  “S-Curve,”  3^(6),  as 

5a  (5)  =  Rl  (s-  I)  Tc  -Rl  (5  +  I)  Tc  (56) 

which  represents  the  loop’s  tracking  curve.  The  expected  loop  discriminator  output  is  consequently 

e{t,6)  ='^K,PS^{6)  (57) 

where  8  is  the  NCDLL  tracking  error. 

Recall  that  Equation  4  is  the  autocorrelation  function  for  a  maximal  length  PRN  sequence. 
Upon  substitution  of  Equation  4  into  Equation  56,  and  after  some  simplification  (10),  the  S-Curve 
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IS  given  as 


0  for-iV  +  l  +  f<5<-(l  +  f) 

-  [1  +  (l  +  ;^)  (6  +  f)]  for  -  (l  +  f)  <  <5  <  “f 

-  [l  -  (l  +  :^)  ((5  +  f)]  for-f  <  5  <  -  (l  -  f) 

^  2  (1  +  i)  [2  -  (1  +  i)  A]  ^  for-  (1  -  f)  <  ^  <  (1  -  f)  (58) 

[1  +  (1  +  7f)  -  t)Y  -  W  for  (1  -  f )  <  ^  <  f 

[1  -  (1  +  jf)  (^  -  f)]^  -  forf  <  5  <  1  +  f 

0  forl+f<«<iV-l-f 

for  A  >  1.0,  or 

0  for-JV  +  l  +  f<«<-(H-f) 

-  [l  +  (^  +  y)  (l  +  ;^)]  ^  for  -  (l  +  y)  <  5  <  (f  -  1) 

-2  (1  +i)  A  [1  +  (1  +  i)  for-(f  -  1)  <  ^  <  -f 

2(1  +  i)  [2-  (l+i)  A]«  (59) 

2  (1  +  i)  A  [1  -  (1+  i)  «|  fotf  <  «  <  (1  -  f) 

[1  -  (1  +  w)  -  f)]^  -  W  for  (1  -  f )  <  ^  <  (1  +  f ) 

0  forl  +  f<5<iV-l-f 

for  A  <  1.0. 

Recall  from  Chapter  I  that  the  autocorrelation  function  of  a  periodic  maximal  length  spread¬ 
ing  code  of  period  N  is  also  periodic  with  period  N.  Therefore,  the  S-Curve  defined  above  is 
periodic  with  period  N.  Equations  58  and  59  give  principal  function  values  which  repeat  every  N 
chips.  Plots  of  Equations  58  and  59  are  shown  in  Figures  14  and  15,  respectively. 

It  is  important  to  note,  as  evident  in  Figures  14  and  15,  that  the  slope  of  the  S~curve  near 
(5  =  0  is  dependent  upon  the  correlator  spacing,  A.  This  slope  is  zero  for  a  correlator  spacing  of 
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S-Curve 


Figure  14  The  Non-Coherent  Delay  Lock  Loop  S-Curve  for  A  =  1  Chip 

2  chips;  therefore,  a  correlator  spacing  of  less  than  2  chips  is  always  used.  Additionally,  note  from 
Figures  14  and  15  that  the  size  of  linear  tracking  region  about  <5  =  0  decreases  with  decreasing  A. 
A  larger  linear  tracking  region  aids  in  code  acquisition  (prior  to  code  lock  and  tracking)  because 
the  tracking  loop  is  able  to  tolerate  larger  tracking  errors  while  establishing  code  lock.  For  this 
reason,  many  receivers  employ  NCDLL’s  with  variable  correlator  spacing.  During  code  acquisition, 
the  correlator  spacing  is  set  at  more  than  1  chip.  Then,  once  code  lock  is  achieved,  the  correlator 
spacing  is  decreased  to  some  smaller  value  (usually  1  chip  for  a  typical  receiver  or  0.1  chip  for  a 
narrow  correlator  receiver). 

3. ^.0.2  Additive  White  Gaussian  Noise  NCDLL  Analysis.  Having  determined  the 
signal  component  of  the  discriminator  output,  next  consider  the  component  due  to  noise.  Recall 
that  the  noise  components  at  the  correlator  outputs  were  given  in  Equations  50  and  51.  Referring 
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S-Curve 


Figure  15  The  Non-Coherent  Delay  Lock  Loop  S-Curve  for  A  =  0.5  Chip 
back  to  those  equations,  define  the  noise  as  follows. 


I  ^  /  A  \ 

(60) 

nL,in(*)  =  -  fd  -  yTc)  n'{t) 

(61) 

where  the  “E”  and  “L”  subscripts  denote  the  early  and  late  correlator  arms,  respectively,  the  “in” 
subscript  denotes  noise  at  the  NCDLL  input,  and 

n\t)  —  \/2ni{t)  cos(a;/i?i  -  <j6')  -  V2nQ{t)  sin  -  0')  (62) 

Next,  to  find  the  noise  component  at  the  bandpass  filter  output,  the  noise  power  spectral  density  at 
the  input  is  needed.  Let  Sn^  and  Sc  be  the  power  spectral  densities  of  the  noise  process,  n'(t),  and 
the  spreading  code,  c{f),  respectively,  and  *  represent  the  convolution  operator.  Since  the  noise 
and  the  spi'eading  code  are  assumed  to  be  independent,  the  power  spectral  densities  of  and 
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nL,in{t)  are  found  by  convolving  the  spreading  code  power  spectral  density  with  that  of  n'{t) 

Snj,u.  (w)  =  ^  Sn'  (o;)  *  Sc  (a>)  j  ^  E,L  (63) 

where  j  designates  either  the  early  or  late  arm  of  the  NCDLL. 

Since  the  spreading  code  has  a  large  bandwidth,  the  effect  of  the  convolution  of  Equation  63 
is  to  spread  the  noise  power  over  a  bandwidth  wider  dhan  the  original.  The  bandpass  filter  has 
a  smaller  bandwidth  with  respect  to  the  spreading  code  and  only  passes  a  small  fraction  of  the 
(spread)  noise  power  spectrum.  Consequently,  the  power  spectrum  in  Equation  63  only  needs  to 
be  evaluated  at  frequencies  near  the  intermediate  frequency,  wjp  (10).  Thus,  Equation  63  can  be 
approximated  by 


(w)  «  (w/f) 


^  (^)  |w  ±  Wif\  <  itBn 
0  otherwise 


(64) 


where  Bn  is  the  one  sided  bandwidth  of  the  noise  process,  Ki/2  accounts  for  the  RF-IF  conversion 
loss  and  the  power  division,  and  No/2  is  the  power  spectral  density  magnitude  of  n  (t),  for  all  co. 

With  the  approximate  noise  power  spectrum  given  in  Equation  64,  the  output  of  the  bandpass 
filters  can  be  calculated  using  the  relationship 


5o,.t(a;)  Siniu;) 


(65) 


where  Sout  (^)  and  Sin{uj)  are  the  power  spectra  at  the  output  and  input,  respectively,  of  a  filter  with 
transfer  function  H{cj).  Assuming  that  the  bandpass  filter  is  an  ideal  filter  with  noise  bandwidth 
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Bn  yields  the  following  power  spectrum  for  the  noise  at  the  filter  output. 


^  |lp  for  |w  ±  W/fl  <  ttSjv 

(66) 

=  “i^iiVo  for  \co  ±  ujif\  <  '^Bn 

With  the  noise  at  the  bandpass  filter  output  thus  characterized,  the  discriminator  output  due  to 
the  signal  and  the  noise  can  be  calculated  as 


e{t,S)  =  [ZsL{t)  +  ZnL{t)fip  -  [ZsEit)  +  ZnE  {t)fip  (67) 

where  Zsl  (i)  and  Zse  (0  signal  components  at  the  output  of  the  late  and  early  arm  bandpass 

filters,  respectively,  ZnL  and  ZnE  are  the  corresponding  noise  components,  and  the  Ip  denotes  the 
lowpass  component  of  the  signals. 

It  can  be  shown  (10)  that  the  discriminator  output  of  Equation  67  is  equivalently 


€{t,6) 


[p-t) 


+  y/2^!^Rl  [(5  -  I)  Te  -  Rl  [(^  +  y)  T, 

X  cos[<^  -  <f>'  +  9d{t  -  7d)] 

+  yf2^i^Rl  [(5  -  I)  Tc]  n'x,,Q(f)  -  Rl  [(5  +  I)  Te] 


(68) 


X  sin[(^  -  0d{t  -  Td)] 


where  rtEi^pi,  and  are  the  in-phase  and  quadrature  components  of  the  noise 

signals  in  the  early  and  late  correlator  arms,  respectively. 
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If  the  noise  components  of  the  discriminator  output  are  denoted  by  the  discriminator 

of  Equation  68  can  be  written  more  succinctly  as 

e{t,8)  -  ^K^PS^{6)  +  n,{t)  (69) 

where  S^{8)  is  the  loop  S-Curve  defined  in  Equation  56  and  6  was  previously  defined  to  be  the 
normalized  tracking  error,  6  =  —  Td)  fTc- 

Assuming  no  data  modulation  is  present,  the  NCDLL  can  be  represented  (10)  by  the  equiv¬ 
alent  model  shown  in  Figure  16.  Although  it  is  unrealistic  to  assume  no  data  modulation  exists, 
this  assumption  leads  to  a  simpler  model.  Additionally,  this  assumption  leads  to  a  worst  case 
performance  model  in  the  sense  that  this  model  will  have  the  maximum  possible  noise  component 
at  the  discriminator  output.  Therefore,  this  model  is  valuable  for  analyzing  the  performance  of  the 
NCDLL. 


(t) 


Figure  16  The  Non-Coherent  Delay  Lock  Loop  Non-Linear  Model 

In  Figure  16,  the  voltage  controlled  oscillator  (VCO)  is  described  by  the  following  input/output 
relationship. 
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where  gc  is  the  VCO  gain  with  units  of  Hz/V.  Noting  from  Figure  16  that  u(A)  is  the  loop  filter 
output  and  that  the  loop  filter  input  is  the  discriminator,  Equation  70  can  be  re-written  as 


MO 

Tc 


“  a)dadX 

Jq  J —oo 


(71) 


where  t;(A)  from  Equation  70  has  simply  been  replaced  by  the  convolution  of  the  discriminator 
output,  and  the  loop  filter  impulse  response,  f{t). 

Substituting  from  Equation  69  yields 


Tc 


=  gc 


ft  [i7iriP5A(^(«))  +  ne(a) 

Jo  J-oo 


/(A  —  a)dad\ 


(72) 


Equation  72  can  be  linearized  if  the  tracking  error,  6^  is  small.  Under  this  assumption,  the 
non-linear  S-curve  can  be  replaced  with  the  following  linear  approximation  (10) 


S^{8) 


(73) 


where  Equation  73  comes  from  Equation  56  evaluated  in  the  region  about  ^  =  0. 

Replacing  the  non-linear  element  5a  (^)  in  Figure  16  with  the  approximation  of  Equation 
73  yields  the  linear  model  of  Figure  17.  In  this  model,  the  S-curve  element,  5a  (^),  and  the  gain 
element,  \KiP^  from  Figure  16  are  replaced  by  a  single  equivalent  gain  element,  where 


Kd 


(74) 
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Figure  17  The  Non-Coherent  Delay  Lock  Loop  Linear  Model 

Let  F(s)  denote  the  Laplace  transform  of  f(t).  Using  the  linear  model  of  Figure  17,  the  closed 
loop  transfer  function,  H(s),  for  the  NCDLL  is 


H{s) 


Kd9cF(s) 
s  +  KdgcF{s) 


Td(s) 


(75) 


Additionally,  the  Laplace  transform  of  the  tracking  error  can  be  shown  (10)  to  be 


6{s) 


rdjs)  r _ s _ ■ 

Tc  [s  +  isrrfffci^(s). 


(76) 


Finally,  the  root-mean-square  (RMS)  tracking  jitter  for  a  loop  with  two-sided  noise  bandwidth  of 
Wl  is  given  by  (10) 


<^6  = 


2^1 


Wl 


(77) 


where  77/2  is  the  discriminator  noise  component  power  spectral  density,  approximately  given  as  (10) 


V 

2 


^{KiNofBN  + 


Tc 


(78) 
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S.  5  Conclusion 


This  chapter  describes  the  complex  cepstrum  adaptive  filter.  It  begins  with  the  input  signal 
model,  which  consists  of  a  direct  path  component,  and  a  delayed  and  attenuated  refiection  of 
the  direct  path.  Next,  each  component  of  the  CCAF  is  described.  Detail  is  given  concerning 
the  function  and  design  of  the  CCAF.  Finally,  an  analysis  of  the  non-coherent  delay  lock  loop  is 
presented.  The  non-linear  and  linear  models  of  the  NCDLL  are  given.  The  chapter  concludes  by 
presenting  the  NCDLL  transfer  function,  tracking  error,  and  RMS  tracking  jitter. 
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IV.  Computer  Simulations  and  Results 


4.1  Overview 

This  chapter  presents  computer  simulations  of  the  theory  presented  in  the  previous  chapters. 
The  purpose  of  these  simulations  is  to  verify  the  functionality  of  the  complex  cepstrum  adaptive 
filter  for  multipath  mitigation,  to  characterize  the  noise  performance  of  the  CCAF,  and  to  compare 
the  CCAF  performance  to  that  of  a  narrow  correlator  non-coherent  delay  lock  loop.  Two  series  of 
simulations  are  run  as  follow: 

1.  CCAF  noise  performance  characterization.  Simulations  are  run  with  constant  multipath 
parameters  and  various  SNR  levels  to  determine  the  noise  sensitivity  of  the  CCAF. 

2.  CCAF  performance  comparison  with  a  narrow  correlator  NCDLL.  These  simulations 
are  run  using  the  CCAF  in  conjunction  with  a  standard  non-coherent  delay  lock  loop, 
employing  a  1.0  chip  correlator  spacing.  The  performance  of  the  CCAF  is  compared 
with  that  of  a  narrow  correlator  NCDLL  employing  a  0.1  chip  correlator  spacing. 

All  simulations  are  written  for  Mat  lab  version  5  and  run  on  the  Sun  workstations  available 
at  the  Air  Force  Institute  of  Technology  or  on  a  stand  alone  personal  computer.  Simulation  results 
are  analyzed  using  Matlab.  All  Matlab  function  files  used  for  the  simulations  are  included  in 
Appendix  A. 

4.2  Simulation  Parameters 

For  all  simulations,  an  attempt  is  made  to  make  all  parameters  consistent  with  actual  GPS 
parameters.  In  cases  where  actual  GPS  receiver  parameters  are  unknown,  simulation  parameters 
are  chosen  to  be  values  which  could  reasonably  be  expected  for  a  GPS  receiver.  Additionally,  when 
some  freedom  of  choice  is  available,  parameters  are  chosen  to  satisfy  certain  desirable  qualities  to 
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enhance  simulations  (for  instance,  decreased  run  times).  The  subsections  that  follow  discuss  the 
simulation  parameters  chosen. 

4.2.1  DS/SS  Spreading  Code.  Actual  GPS  C/A  spreading  code  is  used  for  these  simula¬ 
tions,  The  code  is  generated  using  the  GPS  toolbox  for  Matlab,  available  from  the  Electrical  and 
Computer  Engineering  Department  at  AFIT.  The  code  period  is  1023  chips  long  (approximately  1 
ms).  Each  chip  is  divided  into  100  samples  to  achieve  a  resolution  of  0.01  chips.  Thus,  the  sampling 
frequency  for  all  simulations  is  102.3  MHz. 

4.2.2  Carrier  and  Direct  Path  Signals.  Since  actual  GPS  receivers  typically  operate  at  an 
intermediate  frequency,  the  “carrier”  frequency  for  all  simulations  is  an  appropriate  intermediate 
frequency.  In  order  to  keep  vector  sizes  manageable,  thus  reducing  run  times,  the  intermediate 
frequency  for  these  simulations  is  chosen  to  be  5  times  the  chipping  rate,  or  5.115  MHz.  For 
convenience,  the  carrier  is  chosen  to  have  a  power  of  1/2  W  at  the  receiver.  Also  for  convenience,  the 
carrier  phase,  Oq,  and  direct  path  propagation  delay,  r^,  are  both  chosen  to  equal  zero.  Finally,  the 
direct  path  attenuation,  ao  is  chosen  to  equal  1.  The  choices  of  direct  path  phase,  propagation  delay, 
and  attenuation  factor  have  the  effect  of  normalizing  the  multipath  parameters  to  the  direct  path 
values.  Thus,  these  simulations  could  represent  any  actual  values  for  these  direct  path  parameters 
if  the  multipath  parameters  are  scaled  by  appropriate  values. 

4.2.S  Multipath  Parameters.  As  mentioned  in  the  previous  paragraph,  the  direct  path 
parameters  are  chosen  such  that  the  multipath  signal  will  be  normalized  to  the  direct  path.  As  such, 
the  multipath  attenuation  factor,  ai  is  chosen  to  vary  between  zero  and  one.  For  most  simulations, 
0.2  <  ai  <  0.8  is  used.  Additionally,  the  multipath  delay,  is  varied  between  approximately 
10  nanosecond  and  1.5  microseconds.  These  delays  correspond  to  a  normalized  delay  (in  chips)  of 
0.1  <  a  <  1.5. 
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NCDLL  Design.  The  non-coherent  delay  lock  loop  used  for  these  simulations  is 
an  open-loop  design.  This  open-loop  design  uses  the  NCDLL  S-curve,  discussed  in  Chapter  III,  to 
characterize  the  tracking  bias  present  after  processing  by  the  CCAF.  The  CCAF  tracking  bias,  using 
an  NCDLL  with  standard  1  chip  correlator  spacing  is  compared  with  the  tracking  bias  of  a  narrow 
correlator  NCDLL,  with  0.1  chip  spacing,  which  does  not  employ  CCAF  processing.  The  open  loop 
design  is  chosen  because  tracking  bias  is  the  primary  metric  for  comparing  one  multipath  mitigation 
technique  to  another.  Additionally,  since  this  thesis  focuses  on  a  signal  processing  technique  used 
prior  to  tracking  in  the  NCDLL,  and  not  on  processing  techniques  within  the  NCDLL  itself,  the 
performance  of  the  NCDLL  is  not  in  question  for  this  work.  Rather,  the  amount  of  improvement 
provided  by  the  pre-processing  of  the  CCAF  is  the  parameter  to  be  studied.  Therefore,  closed-loop 
analysis  of  the  NCDLL,  to  study  parameters  such  as  tracking  jitter,  is  not  important  for  this  work. 

4.2.5  Simulation  Parameters  Summary.  The  simulation  parameters  from  the  previous 
subsections  are  summarized  below. 

•  Spreading  Code:  Actual  GPS  C/A  code 

•  Spreading  Code  Period:  1023  chips 

•  Sampling  Frequency:  102.3  MHz  (100  samples/chip) 

•  Carrier  Intermediate  Frequency:  5.115  MHz  (5  times  chipping  rate) 

•  Carrier  Received  Power:  1/2  W 

•  Direct  Path  Phase,  ^o-  0  rad 

•  Direct  Path  Propagation  Delay,  Td'.  0  sec 

•  Direct  Path  Attenuation  Factor,  aQi  1 

•  Multipath  Reflection  Attenuation  Factor,  ai:  0.2  <  o-i  <  0.8 

•  Multipath  Reflection  Normalized  Delay,  a:  0.1  <  o  <  1.5  chips 


•  NCDLL  Correlator  Standard  Spacing:  1  chip 


•  NCDLL  Correlator  Narrow  Spacing:  0.1  chip 

Signal-to-Noise  Ratio  Analysis 

SignaLto-noise  ratio  is  a  fundamental  concern  for  all  communication  systems.  SNRs  must  be 
sufficiently  high  to  allow  receivers  to  perform  the  necessary  processing  to  extract  the  information 
from  the  signal.  For  GPS  C/A  code,  the  typical  SNR  is  -14.9  dB  at  the  spectral  peak  (9).  For 
most  communication  systems,  a  negative  SNR  would  prevent  communication;  however,  for  spread 
spectrum  systems,  such  as  GPS,  the  processing  gain  allows  for  communication  even  when  negative 
SNRs  exist.  In  this  research,  though,  it  is  found  that  cepstral  processing  requires  extremely  high 
SNR  prior  to  any  processing  gain  being  achieved. 


Noisy  and  Ideal  Complex  Cepstra 


Figure  18  Complex  Cepstrum  of  Noisy  and  Ideal  Signal 


For  the  signal  model  presented  in  Chapter  III,  simulation  results  show  that  the  SNR  must  be 
approximately  225  dB  or  higher  to  achieve  good  signal  reconstruction  for  the  direct  path  estimate. 
Figure  18  demonstrates  the  effects  of  AWGN  on  the  complex  cepstrum.  In  Figure  18  the  complex 
cepstrum  of  a  signal  in  the  absence  of  AWGN  is  compared  to  that  for  a  signal  in  AWGN  with 
an  SNR  of  100  dB.  It  is  seen  that  the  noise  causes  the  complex  cepstrum  to  appear  very  noisy 
even  with  such  a  high  input  SNR.  This  random  noise  in  the  complex  cepstrum  causes  errors  in  the 
multipath  filtering  process  which  lead  to  errors  in  the  recovered  signal.  Due  to  the  randomness  of 
the  data  points  in  the  complex  cepstrum,  the  filtered  values,  which  are  the  average  of  the  adjacent 
points,  also  have  a  certain  randomness.  This  has  the  effect  of  causing  multipath-like  amplitude 
changes  in  the  recovered  signal.  Plots  demonstrating  this  effect  are  shown  in  Figures  19  and  20 
.  These  amplitude  changes  greatly  degrade  or  even  prevent  tracking  of  the  received  signal  by  the 
code  tracking  loop.  For  this  reason,  the  complex  cepstrum  is  found  to  be  applicable  for  code  loop 
multipath  mitigation  only  under  unreasonably  high  SNR  conditions. 
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Figure  19  Recovered  Signal  from  a  Noisy  Complex  Cepstrum  Process 
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Recovered  Signal 


Figure  20  Recovered  Signal  from  a  Noisy  Complex  Cepstrum  Process  (Expanded  Time  Scale) 

The  SNR  threshold  effect  is  demonstrated  in  Figure  21.  In  this  figure,  the  mean  squared 
amplitude  error  of  the  recovered  signal  is  plotted  on  a  semilog  scale  against  the  input  SNR.  Figure  21 
shows  an  exponential  increase  in  mean  squared  amplitude  error  for  SNRs  which  are  below  160  dB.  In 
this  region,  the  recovered  signal  is  severely  degraded,  preventing  tracking  by  the  NCDLL.  However, 
for  good  tracking  performance  the  SNR  threshold  is  somewhere  near  225  dB,  considerably  higher 
than  the  160  dB  threshold.  This  is  because  even  very  small  differences  between  the  recovered  signal 
and  the  actual  direct  path  signal  lead  to  poor  tracking  performance  by  the  NCDLL,  as  compared 
to  a  NCDLL  tracking  the  multipath  signal.  Figure  21  shows  that  the  mean  squared  amplitude 
error  begins  to  increase  for  SNR  levels  below  approximately  225  dB.  Although  the  MSE  is  still 
fairly  small  in  the  region  around  225  dB  SNR,  it  is  large  enough  that  the  cepstrum  processing  leads 
to  tracking  biases  which  are  larger  than  those  for  the  same  non-coherent  delay  lock  loop  tracking 
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the  multipath  signal  without  cepstrum  processing.  For  this  reason,  the  input  SNR  must  be  above 
approximately  225  dB  to  achieve  good  multipath  mitigation  using  cepstrum  processing. 


Figure  21  SNR  Threshold  Effect 


To  improve  noise  performance,  an  input  filter  would  typically  be  used.  However,  in  this  case, 
filtering  the  input  signal  causes  severe  degradation  of  the  recovered  signal.  To  understand  why 
this  is  so,  refer  back  to  the  basic  signal  model  of  Equations  29  and  30.  Recall  that,  using  these 
signal  models,  and  following  the  derivation  for  calculating  the  complex  and  power  cepstra,  these 
cepstra  are  actually  time  domain  deconvolution  tools.  If  the  signal  is  filtered  prior  to  calculating 
the  cepstra,  the  model  of  Ec|uations  29  and  30  no  longer  applies.  The  new  model  must  account  for 
the  fact  that  the  signal  has  been  convolved  with  the  time  domain  impulse  response  of  the  filter. 
This  second  convolution  (filtering),  changes  the  complex  cepstrum,  and  prevents  proper  multipath 
mitigation  using  the  technicpies  described  in  Chapters  II  and  III.  Figure  22  shows  a  plot  of  the 
complex  cepstrum  of  a  multipath  signal  calculated  before  and  after  bandpass  filtering  the  time 
domain  signal  to  demonstrate  this  point.  The  signal  in  Figure  22  has  no  noise  added  so  that  the 
effect  is  due  solely  to  the  input  filter. 
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Typical  Complex  Cepstrum 


Figure  22  Bandpass  Filtering  Effect  on  the  Complex  Cepstrum 


Cepstrum  Processing  Analysis 

The  cepstrum  processing  for  multipath  mitigation  consists  of  two  steps,  as  described  in  Chap¬ 
ter  III.  First,  the  power  cepstrum  is  used  to  detect  the  multipath  signal.  Then,  the  complex  cep¬ 
strum  is  used  to  filter  out  the  multipath.  Although  both  these  steps  are  sensitive  to  the  input 
SNR,  the  minimum  SNR  is  set  using  the  complex  cepstrum.  The  power  cepstrum  is  significantly 
less  sensitive  to  noise  effects  than  the  complex  cepstrum.  This  observation  will  be  discussed  more 
in  the  section  that  follows. 

Multipath  Detector.  As  described  in  Chapter  III,  the  power  cepstrum  magnitude 
can  be  used  to  detect  the  multipath  delay  by  searching  for  the  first  delta  funtion  in  the  power 
cepstrum.  Simulation  shows  this  process  is  relatively  immune  to  noise  effects  assuming  the  input 
SNR  is  about  30  dB  or  higher,  and  the  threshold  for  declaring  a  multipath  peak  is  set  appropriately. 
Magnitude  plots  of  the  power  cepstrum  for  typical  GPS  multipath  signals  are  shown  in  Figures  23 
and  24. 
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Typical  Power  Cepstrum  Magnitude  Plot 


Figure  23  Power  Cepstrum  for  a  GPS  Multipath  Signal,  SNR  =  30  dB,  Multipath  Amplitude 
=  0.5,  Multipath  Delay  =  0.5 

Analysis  of  Figures  23  and  24  reveals  two  important  observations.  First,  the  power  cepstrum 
of  a  GPS  multipath  signal  contains  “rolling  hill”  type  peaks  at  low  quefrencies.  These  rolling  hill 
peaks,  although  not  caused  by  multipath,  sometimes  cause  the  multipath  delay  to  be  incorrectly 
detected.  This  effect  will  be  observed  later.  The  second  observation  is  that,  as  the  SNR  decreases, 
the  magnitude  of  the  delta  functions  in  the  power  cepstrum  decreases.  This  decrease  drives  the 
detection  threshold.  As  the  delta  function  peaks  become  smaller,  it  becomes  more  difficult  to 
accurately  detect  the  multipath  peaks.  For  the  simple  detector  used  in  this  work,  the  detection 
SNR  threshold  is  approximately  30  dB  SNR. 

To  demonstrate  that  cepstral  processing  can  be  used  for  multipath  mitigation,  a  series  of 
simulations  tests  the  cepstrum  system  of  Chapter  III  over  a  variety  of  multipath  scenarios,  for 
SNRs  which  are  above  the  threshold  mentioned  earlier  in  this  chapter.  For  these  simulations,  the 
multipath  delay  is  varied  from  0.1  chip  to  1.5  chips,  for  multipath  amplitudes  equal  to  0.2,  0.4,  0.6, 
and  0,8,  relative  to  the  direct  path  amplitude.  For  each  pair  of  multipath  delays  and  amplitudes. 
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Typical  Power  Cepstrum  Magnitude  Plot 


Figure  24  Power  Cepstrum  for  a  GPS  Multipath  Signal,  SNR  ==  250  dB,  Multipath  Amplitude 
=  0.5,  Multipath  Delay  =  0.5 

10  trials  are  run  adding  independent  white  Gaussian  noise  to  the  input  signal.  Plots  of  the  detected 
multipath  delay  (dots)  versus  the  actual  multipath  delay  (line)  are  shown  in  Figures  25  through  36. 
Four  observations  should  be  made  concerning  these  plots.  First,  the  detector  works  well  over  a  wide 
variety  of  multipath  scenarios.  Three  points  are  observed  when  the  multipath  detector  performance 
is  not  ideal,  corresponding  to  delays  of  0.1  chip,  0.3  chip,  and  1.0  chip.  The  first  two  problem  spots 
are  due  to  the  rolling  hill  peaks  in  the  power  cepstrum  at  these  low  ciuefrencies.  The  problem  at 
1.0  chip  is  actually  the  way  the  detector  is  designed.  As  discussed  in  Chapter  III,  the  power  and 
complex  cepstrum  both  contain  natural  peaks  at  the  chip  interval.  Simulations  demonstrated  that 
these  peaks  are  important  for  signal  reconstruction,  and  should  not  be  filtered  out.  Therefore, 
since  the  adaptive  weighted  comb  filter  is  designed  not  to  filter  at  the  chip  interval,  the  multipath 
detector  is  designed  not  to  detect  multipath  delays  of  1,0  chip.  This  design  has  the  weakness  that 
when  a  multipath  reflection  exists  at  exactly  1.0  chip  delay,  this  cepstrum  process  cannot  remove 
the  multipath.  This  weakness  will  be  demonstrated  in  the  next  section.  Finally,  when  the  multipath 
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reflection  delay  is  incorrectly  detected,  the  detected  delay  is  always  2.0  chips.  This  is  due  to  the 
first  harmonic  of  the  natural  peak  at  1.0  chip  delay  being  detected  as  the  multipath  signal,  and  is 
due  to  the  fact  that  the  detector  is  designed  with  the  assumption  that  a  multipath  reflection  exists 
at  some  delay. 


Detected  vs.  Actual  Multipath  Delay 


Figure  25  Multipath  Detector  Test  Results  for  SNR  =  1000  dB  and  Multipath  Amplitude  =  0.2 
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Detected  vs.  Actual  Multipath  Delay 


ipath  Detector  Test  Results  for  SNR  250  dB  and  Multipath  Amplitude  =  0.8 


Detected  vs.  Actual  Multipath  Delay 


[path  Detector  Test  Results  for  SNR  =  223  dB  and  Multipath  Amplitude  =  0.2 
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^.^.2  Multipath  Mitigation.  Following  multipath  detection,  the  complex  cepstruin  and  the 
adaptive  weighted  comb  filter  are  used  to  remove  the  reflected  signal  in  the  cepstral  domain.  Then, 
a  direct  path  estimate  is  obtained  through  an  inverse  complex  cepstrum  process.  To  determine  the 
effectiveness  of  the  complex  cepstrum  filtering  process,  two  figures  of  merit  are  used.  First,  the 
mean  squared  amplitude  error  between  the  direct  path  estimate,  and  the  actual  direct  path  signal 
is  calculated  from  the  following  formula 

MSE  =  i  ('^9) 

n  =  0 

where  x(n)  and  x{n)  are  the  data  points  of  the  direct  path  signal  and  the  direct  path  estimate, 
respectively. 

The  second  figure  of  merit  is  the  non-coherent  delay  lock  loop  tracking  bias.  This  bias 
is  calculated  from  the  S-curve  presented  in  Chapter  III.  The  tracking  bias  is  defined  to  be  the 
magnitude  of  the  delay  or  advance  of  the  recovered  signal  S-curve  relative  to  the  direct  path  S- 
curve.  The  advance  or  delay  is  measured  at  the  point  where  the  S-curves  are  equal  to  zero.  This  is 
illustrated  in  Figure  37  which  shows  a  typical,  noise-free  S-curve  plot  with  the  appropriate  tracking 
biases  noted  on  the  figure. 

Mean  Squared  Error  Analysis.  Figures  38  through  41  show  the  average 
mean  squared  amplitude  eri*or  in  the  direct  path  estimate  for  multipath  amplitudes  of  0.2,  0.4,  0.6, 
and  0.8,  and  multipath  delays  from  0.1  chip  to  1.5  chips,  both  relative  to  the  direct  path  signal. 
For  each  pair  of  multipath  amplitudes  and  delays,  10  trials  are  run  at  each  of  three  different  SNR 
levels,  1000  dB,  250  dB,  223  dB;  the  figures  show  the  average  values.  The  SNRs  are  chosen  to 
collect  data  well  above,  slightly  above,  and  near  the  threshold  previously  discussed.  This  gives  an 
idea  of  how  sensitive  the  complex  cepstrum  filtering  process  is  to  input  SNR,  assuming  the  SNR  is 
above  the  threshold. 
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Mean  Squared  Amplitude  Errror  of  Recover 


Figure  39 


Figure  40 


Direct  Path  Estimate  Mean  Squared  Error  For  Multipath  Ampli 


Mean  Squared  Amplitude  Errror  of  Recovered  Signal 


Direct  Path  Estimate  Mean  Squared  Error  For  Multipath  Ampli 
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Figure  41  Direct  Path  Estimate  Mean  Squared  Error  For  Multipath  Amplitude  —  0.8 

Examination  of  Figures  38  through  41  reveals  that  the  mean  squared  error  is  reasonably 
insensitive  to  SNR,  assuming  the  threshold  condition  is  met.  In  Figure  41,  a  condition  where  the 
SNR  threshold  is  not  met  is  observed.  In  this  plot,  the  amplitude  MSE  for  a  multipath  signal 
delayed  0.1  chip  relative  to  the  direct  path,  with  an  input  SNR  of  223  dB  is  extremely  large.  This 
is  a  direct  result  of  the  SNR  threshold  not  being  met  for  that  pair  of  multipath  parameters.  In 
the  next  subsection,  it  will  be  seen  that  this  large  mean  squared  amplitude  error  in  the  estimated 
direct  path  signal  amplitude  leads  to  a  large  tracking  bias.  This  figure  also  shows  that  223  dB  is 
near  the  SNR  threshold  also  for  a  multipath  delay  of  0.15  chips.  Other  than  those  two  data  points, 
the  mean  squared  error  performance  is  adequate  for  all  other  delays  at  all  three  SNRs. 

^.4‘2.2  TracMng  Bias  Analysis.  The  second  figure  of  merit  for  characterizing  the 
cepstrum  multipath  mitigation  process  is  the  tracking  bias.  To  perform  this  analysis,  simulations 
are  run  for  multipath  amplitudes  of  0.2,  0.4,  0.6.  and  0.8  at  SNRs  of  1000  dB,  250  dB,  and  223  dB. 
For  each  multipath  amplitude  and  SNR,  the  multipath  delay  is  varied  from  0.1  to  1.5  chips.  Again, 
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10  trials  are  run  for  each  set  of  multipath  parameters,  and  the  data  presented  is  the  average  for 
those  10  experiments.  To  characterize  the  cepstrum  performance,  the  tracking  bias  of  the  cepstrum 
process  is  calculated  and  compared  to  the  tracking  bias  for  a  narrow  correlator  non-coherent  delay 
lock  loop  employing  0.1  chip  correlator  spacing.  This  comparison  is  performed  to  characterize  the 
cepstrum  process  in  relation  to  a  cuiTent,  commercially  available,  method  of  multipath  mitigation. 
The  performance  of  the  cepstrum  process  is  compared  to  that  of  the  narrow  correlator  in  Figures 
42  through  53. 


Figure  42  Tracking  Bias  Comparison  for  SNR  =  1000  dB  and  Multipath  Amplitude  =  0.2 
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Multipath  Delay 


Bias  Comparison  for  SNR  = 


1000  dB  and 


Tracking  Bias  vs.  Multipath  Del. 


Examination  of  Figures  42  through  53  reveals  that  the  cepstrum  process  generally  exhibits 
lower  tracking  bias  than  the  narrow  correlator,  when  the  SNR  threshold  requirement  is  met.  Ex¬ 
ceptions  to  this  general  rule  occur  when  the  multipath  delay  is  small  and  when  it  equals  1  chip. 
For  small  multipath  delay,  the  adaptive  comb  filter  must  filter  the  complex  cepstrum  at  more  data 
points,  introducing  more  opportunity  for  degradation  in  the  recovered  signal.  For  instance,  if  the 
multipath  delay  is  0.1  chip,  the  complex  cepstrum  must  be  filtered  at  every  integer  multiple  of  0.1 
chip.  If,  on  the  other  hand  the  delay  is  0.5  chip,  it  will  be  filtered  at  every  integer  multiple  of  0.5 
chip.  The  first  example  will  be  filtered  at  5  times  as  many  quefrencies  as  the  second,  introducing 
more  opportunities  for  filtering  error  in  the  recovered  signal.  Additionally,  it  has  been  seen  that, 
at  small  multipath  delays,  the  detector  frequently  misdetects  the  multipath  delay.  This  leads  to 
the  recovered  direct  path  estimate  more  closely  resembling  the  composite  multipath  signal  rather 
than  the  direct  path  signal.  However,  all  current  multipath  mitigation  techniques  exhibit  poor 
performance  for  small  reflection  delays. 

As  previously  mentioned,  the  second  exception  occurs  when  the  reflection  delay  is  1.0  chip. 
Due  to  the  importance  of  the  cepstral  peaks  at  each  integer  multiple  of  1.0  chip,  as  discussed 
previously,  the  adaptive  weighted  comb  filter  is  designed  not  to  filter  the  cepstrum  at  multiples 
of  the  chip  interval.  Therefore,  by  design,  if  the  multipath  signal  is  delayed  exactly  1.0  chip, 
relative  to  the  direct  path,  the  cepstrum  process  performs  no  filtering,  and  the  original  composite 
multipath  signal  is  sent  to  the  tracking  loop.  This  weakness  degrades  performance,  compared  to 
the  narrow  correlator,  in  the  one  case  where  the  multipath  delay  is  1.0  chip;  however,  it  ensures 
better  performance  at  all  other  multipath  delays. 

Finally,  examination  of  Figure  53  reveals  that  the  maximum  tracking  bias  for  the  cepstrum 
process,  when  the  input  SNR  is  223  dB  and  the  multipath  amplitude  is  0.8,  occurs  at  0.1  chip 
multipath  delay.  Recall  that  the  mean  squared  error  of  the  recovered  signal  in  that  case  is  large 
due  to  the  threshold  effect.  The  poor  quality  of  the  recovered  signal  at  this  point  leads  to  poor 
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Table  4  Average  Data  Summary  For  223  dB  SNR 


Parameter 

Multipath  Amplitude 

0.2 

0.4 

0.6 

0.8 

Mean  Squared  Error 

0.0105 

0.0143 

0.0229 

0.0577  (edited) 

Cepstrum  Bias  (chips) 

0.0165 

0.0195 

0.0286 

0.0556 

Narrow  Correlator  Bias  (chips) 

0.0379 

0.0445 

0.0584 

0.1540 

Standard  Correlator  Bias  (chips) 

0.0322 

0.0713 

0.1183 

0.2648 

well  over  a  wide  variety  of  multipath  scenarios  using  very  high  input  SNRs.  The  “edited”  line  in 
Table  4  presents  the  average  mean  squared  error  with  the  one  outlying  data  point  (shown  in  Figure 
41)  removed. 
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V,  Conclusion 


This  thesis  seeks  to  provide  a  new  method  for  GPS  multipath  mitigation.  It  proposes  use  of  the 
cepstrum  to  detect  and  i^emove  multipath  effects  prior  to  code  tracking  by  a  non-coherent  delay 
lock  loop.  In  Chapter  II,  the  theory  of  the  cepstrum  is  presented.  Chapter  III  discusses  the  signal 
model,  the  cepstrum  process,  and  the  non-coherent  delay  lock  loop.  Chapter  IV  presents  the  results 
of  simulations  using  the  cepstrum  for  multipath  mitigation. 

Recall  from  Chapter  I  that  the  objectives  of  this  thesis  are: 

1.  Develop  a  complex  cepstrum  filtering  technique  for  GPS  multipath  interference  removal. 

2.  Modify  a  standard  delay  lock  loop  by  adding  the  filtering  to  the  loop  input. 

3.  Characterize  the  multipath  induced  tracking  error  for  the  modified  delay  lock  loop  in  a 
noiseless  environment. 

4.  Characterize  the  multipath  induced  tracking  error  for  the  modified  delay  lock  loop 
operating  with  typical  GPS  SNRs. 

5.  Compare  and  contrast  the  complex  cepstrum  receiver  performance  to  that  of  a  narrow 
correlator  receiver. 

These  objectives  were  met  through  the  theory  of  Chapters  II  and  III,  and  demonstrated 
through  the  simulations  of  Chapter  IV.  Based  on  the  results  of  Chapter  IV ,  several  conclusions  can 
be  made. 

1.  The  power  cepstrum  can  be  used  to  detect  multipath  interference  prior  to  GPS  code 
tracking. 

2.  The  complex  cepstrum  can  be  used  to  remove  multipath  interference  prior  to  GPS  code 
tracking. 
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3.  Both  the  complex  and  power  cepstra  are  sensitive  to  input  SNR,  with  the  complex 
cepstriim  being  extremely  sensitive  to  input  noise. 

4.  Time  domain  filtering  to  improve  input  SNR  causes  serious  degradations  in  the  cepstral 
domain. 

5.  Due  to  the  noise  sensitivity  of  the  complex  cepstrum,  this  technique  is  not  viable  at  this 
time  for  use  in  actual  GPS  receivers  where  the  input  SNR  is  typically  -14.9  dB. 

Although  the  noise  sensitivity  of  the  complex  cepstrum  makes  it  unusable  for  GPS  multipath 
mitigation,  this  research  shows  that  the  complex  cepstrum  is  a  viable  technique  for  multipath 
mitigation  in  other  circumstances.  This  technique  may  be  applicable  to  other  communication 
systems  with  different  input  signal  structures  and  higher  SNRs.  One  example  of  such  a  system  may 
be  overseas  telephony  where  users  often  hear  echoes  while  talking. 

5.1  Recommendations  for  Future  Research 

Despite  the  fact  that  these  cepstrum  techniques  do  not  appear  viable  for  GPS  multipath 
mitigation  at  this  time,  further  research  is  warranted  to  investigate  possible  improvements  to  make 
the  cepstrum  a  viable  tool  for  GPS  applications.  This  research  includes  the  following: 

1.  Investigate  cepstral  processing  for  multipath  mitigation  within  or  following  a  non¬ 
coherent  or,  perhaps  more  promisingly,  a  coherent  delay  lock  loop,  rather  than  strictly 
prior  to  a  NCDLL.  Correlation  prior  to  cepstral  processing  may  lead  to  the  signal  to 
noise  improvements  necessary  to  make  this  technique  applicable  for  GPS  signals.  A  post 
delay  lock  loop  cepstral  algorithm  has  the  advantage  of  greatly  enhanced  SNR  at  the 
expense  of  a  distorted  measured  autocorrelation  function  (via  a  bank  of  evenly  spaced 
correlators  as  in  the  MEDLL  or  MRDLL  designs).  Cepstral  techniques  could  then  be 
used  to  deconvolve  the  primary  path  autocorrelation  function  from  the  contribution  of 
any  component  reflections. 
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2.  Develop  a  new  multipath  signal  model  that  accounts  for  the  time  domain  convolution 
of  an  input  filter  impulse  response.  Characterize  the  filter  influence  in  the  cepstral 
domain,  and  develop  methods  of  detecting  and  removing  multipath  despite  the  effects 
of  the  input  filter. 

3.  Develop  methods  of  cepstral  domain  filtering  to  minimize  input  noise  effects  prior  to 
multipath  detection  and  mitigation. 

4.  Using  the  theory  of  Chapter  II,  extend  this  research  to  include  multiple  reflections. 

5.  Improve  the  signal  model  to  include  effects  of  the  carrier  tracking  loop,  and  characterize 
the  performance  of  the  cepstrum  process  on  this  signal. 

6.  Characterize  the  cepstrum  process  when  doppler  effects  are  taken  into  account. 

7.  Characterize  the  cepstrum  process  for  GPS  carrier  tracking  multipath  mitigation. 

5,2  Final  Conclusion 

This  work  shows  the  cepstrum  is  viable  for  multipath  mitigation  under  ideal,  noiseless  con¬ 
ditions.  The  simulations  show  that  the  cepstrum  process  requires  an  extremely  high  input  SNR. 
Further  research  is  necessary  to  improve  noise  performance  before  this  will  be  a  viable  technique 
for  GPS  code  tracking  multipath  mitigation. 
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Appendix  A.  Matlab  Function  Files 


A,1  Overview 

The  Appendix  contains  copies  of  all  Matlab  (version  5)  function  filesised  for  simulation  of 
multipath  mitigation  via  cepstral  techniques.  Each  file  is  a  separate  sectkof  the  Appendix.  A 
short  introduction  to  each  section  is  given,  followed  by  a  listing  of  the  fundin  file. 


A. 2  Thesis^SiruMlator  Function  File 

The  function  file  thesis.simulator  is  the  primary  file  used  for  simulatk»f  cepstral  processing 
for  multipath  mitigation.  This  file  creates  the  composite  multipath  signal,  dfeilates  the  power  and 
complex  cepstra,  calls  the  detector  subroutine  to  detect  the  reflection  delay,amoves  the  multipath 
effects  in  the  cepstral  domain,  and  produces  the  direct  path  estimate  thrcsiji  an  inverse  complex 
cepstrum  operation.  Additionally,  this  file  calculates  the  mean  squared  amjAde  error  in  the  direct 
path  estimate,  compared  to  the  actual  direct  path  signal.  Finally,  it  caldttes  the  S-curves  for 
the  direct  path  signal,  the  direct  path  estimate,  the  narrow  correlator  NCBL,  and  the  standard 
correlator  NCDLL.  These  S-curves  are  then  used  to  calculate  the  trackingfeses  of  the  cepstrum 
process,  the  narrow  correlator  NCDLL,  and  the  standard  correlator  NCDOThe  Matlab  function 
file  follows: 

function  [results] =thesis_simulator (am, tm, SNR) 

Vo  THESIS_SIMULATOR  uses  the  cepstrum  to  detect  and  remove  a  nrallpath  reflection  from  a 
7o  composite  signal .  The  direct  path  signal  is  modeled  as  a  co^e  modulated  by  a  BPSK 
7  spreading  code.  The  multipath  reflection  is  modeled  as  an  aajitude  scaled  and  time 
7o  delayed  version  of  the  direct  path  signal.  The  composite  sig^l  is  the  sum  of  the 
%  direct  path  and  reflection  signals. 

7o 

7o  Usage: 

7o  results=thesis_simulator(am,tm,SNR) 

7o 

7o  Where : 

7o  am  =  multipath  reflection  normalized  amplitude 

7p  tm  =  multipath  reflection  normalized  delay  (in  chips) 

7e  SNR  =  the  input  signal  to  noise  ratio  in  dB 
7p 
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7«  This  file  calls  the  following  subroutines: 

7o 

t  expand,  delay,  fpeak,  meanfilt,  scurve,  trackpoint 

I 

7o  This  function  and  all  subroutines  were  written  by  Chuck  Ormsby 

% 

load  cacode  7o  Load  pre-’Stored  C/A  code  generated  using  the  GPS  toolbox  for  Matlab 
fs=100;  %  fs  is  the  sampling  frequency  (ie  100  samples/chip) 
ca=[expand(ca,fs)] ;  %  Produces  100  samples  per  chip 

ca=[ca  zeros(size(ca))]  ;  7o  Zero  pad  the  signal  for  cepstral  processing 
t=linspace(0,2*pi*1023*2,length(ca)) ;  7o  Time  scale  for  the  carrier 
carrier=cos(5*t) ; 

y=ca.*carrier ;  7o  Direct  path  signal 
n=(0:length(ca)-l)/100; 

7®  Add  a  multipath  signal 

d=round(tm*f s) ;  7®  Reflection  delay  Note:  Round  is  necessary  to  make  sure  d  is  an  integer. 
a==am;  7®  Reflection  Amplitude 

mp=a*delay(y  ,d) ;  7®  Multipath  reflection  signal 
cs=y+mp;  7®  Composite  signal 

7®  Convert  the  SNR  in  dB  to  a  noise  amplitude 

na=exp(-SNR/20); 

7®  Generate  a  Gaussian  random  sequence 

randnC ’ state \ sum(  100*clock) )  ;  7®  Reset  the  random  number  generator. 
rv=na*randn(size(cs) ) ; 

7#  Add  the  noise  to  the  signal 
cs=cs+rv;  7®  Noisy  composite  signal 

7®  Caculate  the  power  and  complex  cepstra 

[CS,nd]=cceps(cs) ;  7®  Complex  cepstrum;  nd  is  delay  which  must  be  removed  later 
RCS=rceps(cs)  ;  7®  Power  cepstrum 

7®  Use  the  power  cepstrum  to  detect  the  multipath  delay 
first_peak=fpeak(abs(RCS) ,  .6)  ;  7®  Detect  the  multipath  peak 
detected_delay=(f irst_peak~l)/fs ;  7®  Detected  delay  in  chips 

7®  Filter  the  multipath  out 

max_multiplier=f loor (length(CS)/f irst^peak) ; 
multipliers=l :max_multiplier ; 
f ilter_points=f irst_peak*multipliers ; 

for  index=l : length(f ilter^points) 

if  rem(f ilter_points(index)-*index+l ,  100) "'=1  7®  Do  not  filter  at  the  chip  times 

CS=meanfilt(CS,filter^points(index)-index+l,2) ;  7®  Filters  out  the  delta  functions 
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end 


end 

%  Recover  the  direct  path  estimate 
cal=icceps(CS,nd) ;  %  Inverse  cepstrnm  calculation 

%  Calculate  the  mean  squared  amplitude  error  in  the  recovered  signal 
mse=mean((y-cal)  .  "2)  ;  */«  MSE  in  the  recovered  signal 

7#  Calculate  the  s-curves  of  the  actual  direct  path,  the  direct  path 
%  estimate,  the  narrow  correlator  receiver  and  the  standard  NCDLL 

[sdp,t]=scurve(y,y,fs,l) ;  7o  Direct  path  s-curve  and  tracking  error  scale  (t) 
sc=scurve(cal,y,fs,l) ;  7#  Cepstrum  process  S“curve 
snc=scurve(cs,y,fs,  .1) ;  7«  Marrow  correlator  s-curve 
sncdll=scurve(cs ,y ,f  s,  1)  ;  7o  Standard  correlator  s-curve 

7e  Calculate  the  tracking  biases  for  the  direct  path  estimate,  the 
7p  narrow 

7fl  correlator  receiver,  and  the  standard  NCDLL  relative  to  the  direct  path  signal 

tpdp=trackpoint (sdp) ;  %  Find  the  zero  crossing  point  of  the  direct  path  s-curve 

tpcep=trackpoint (sc) ;  %  Find  the  zero  crossing  point  of  the  cepstriim  process  s-curve 

tpnc=trackpoint (snc) ;  %  Find  the  zero  crossing  point  of  the  narrow  correlator  s-curve 

tpncdll=trackpoint(sncdll);  7o  Find  the  zero  crossing  point  of  the  standard  correlator  s-curve 

cep3ias=(tpcep-tpdp)/fs;  7o  Calculate  the  cepstrum  tracking  bias  in  chips 
nar_bias=(tpnc-tpdp)/fs;  %  Calculate  the  narrow  correlator  tracking  bias  in  chips 
stan_bias=(tpncdll-tpdp)/fs;  7o  Calculate  the  standard  correlator  tracking  bias  in  chips 

7o  Report  the  results  of  the  simulation. 

results= [SNR  am  tm  detected_delay  mse  cep_bias  nar_bias  stan^bias] ; 


A. 3  Expand  Function  File 

Expand  is  the  first  function  file  called  by  thesis  ^Simula, tor.  Expand  expands  an  input  vector 
containing  one  sample  per  time  unit,  to  one  containing  fs  samples  per  time  unit,  where  fs  is 
the  sampling  frequency.  For  this  specific  application,  expand  takes  an  input  vector  of  GPS  C/A 
spreading  code  containing  one  sample  per  chip,  and  expands  the  vector  to  contain  100  samples  per 
chip.  The  Matlab  function  file  follows: 
function  yp=expand2(y, expan) 

%  EXPAND  will  expand  a  vector  by  an  input  factor.  This  expansion  is  equivalent  to  producing 
7o  more  samples  per  time  unit. 
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% 

t  ex.  y=[l  2  3  4] 

7o  yp=expand(y  ,2) 

%  yp=[l  1  2  2  3  3  4  4] 

I 

7o  yp=expand(y, expansion  factor) 

7o 

7o  Written  by  Chuck  Ormsby 
7o  April  20,  1997 
7o 


yp=zeros(l,length(y)*expan) ; 
count=l ; 
count2=l ; 

while  count2<=length(y) 

yp( count : count+expan-l)=y(count2)*ones(l , expan) ; 
count 2= count 2+1 ; 
count=count+expan ; 

end 


A.\  Delay  Function  File 


The  function  file  delay  adds  a  delay  to  the  input  vector  by  prepending  zeros  to  the  beginning 
of  the  vector.  This  function  is  used  to  delay  the  direct  path  signal  vector  for  creation  of  the  reflection 
signal  vector.  The  Matlab  function  file  follows: 

function  y=delay(x,s) 

7#  DELAY  prepends  zeros  to  the  beginning  of  the  input  vector.  This  can  be  thought  of  as 
7p  adding  a  time  delay  to  the  input  vector. 

7o 

7o  y=delay(x,s) 

7c 

7o  y  is  the  delayed  vector,  x  is  the  input  vector,  s  is  the  length  of  the  delay  and  must 
7e  be  positive 
7c 


left=2eros(l ,s) ; 
right=x(l : length (x)-s) ; 
y=  [left  right] ; 
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A, 5  Fpeak  Function  File 


The  function  file  fpeak  detects  the  multipath  reflection  delay.  This  function  compares  each 
data  point  in  the  input  vector  to  the  two  adjacent  data  points  as  described  in  Chapter  III.  When 
a  delta  function  peak  is  detected,  the  function  reports  the  associated  delay  as  the  reflection  delay. 

The  Matlab  function  file  follows: 
function  first_peak=fpeak(x, multiple) 

%  FPEAK  detects  the  first  peak  in  a  function.  It  is  designed  to  work  with  the  CCAF. 
ft  A  peak  is  declared  if  a  point  exceeds  the  points  on  either  side  of  it  by  more  than 

7o  an  input  multiple  which  is  a  percentage  (ie  multiple  =  .5  means  a  point  must  be 

7o  2  times  larger  than  the  adjacent  point).  By  default,  the  first  point  cannot  be  a  peak. 
7o 

7o  first_peak  =  f peak (x, multiple) 

7o 

79  Written  By: 

7o  Chuck  Ormsby 
7«  27  Oct  97 
79 


stop=0; 
index=2; 
while  stop==0 

if  x(index)>=(l/multiple)*x(index-l)  &  x ( index )>=(! /mult iple)*x(index+l) 
if  index"'=101  79  Do  not  declare  an  mp  delay  of  1  chip 
f ir  St  _p  e ak= inde x ; 
stop=l ; 

else 

index=index+l ; 

end 

else 

index=index+l ; 

end 

end 


A. 6  MeanfiU  Function  File 

The  m,eanfilt  function  file  accomplishes  the  removal  of  the  multipath  reflection  effects  from 
the  complex  cepstrum.  This  function  replaces  the  appi'opriate  data  points  in  the  complex  cepstrum 
with  the  average  of  the  adjacent  data  points.  The  Matlab  function  file  follows: 

function  xf=meanfilt (x,m,N) 
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%  MEANFILT  is  a  "mean"  comb  filter  designed  for  use  with  the  complex  cepstrum.  This 
7,  filter  replaces  element  m  of  the  vector,  x,  with  the  mean  of  the  N  surrounding  points. 
7o  N  must  be  even. 

7o 

I  Ex:  x=  [5  2  3  4  5  6] 

7o 

7o  xf=meanfilt(x,3,4) 

7o  x(3)=(5+2+4+5)/4 
I  x(3)=  4 
7,  xf=[5  2  4  4  5  6] 

7o 

7#  Note:  X  and  m  must  be  specified.  If  N  is  not  entered,  a  default  value  of  N=2  is  used. 
7o 

7#  Written  by  Chuck  Ormsby 

7o 

if  nargin==2 
N=2; 

end 

7,  Test  N 

while  rem(N,2)"=0 

N=input(’N  must  be  even.  Please  enter  a  new  N.  O; 

end 

x(m)  =  (sum(x(m-N/2  :m“l)  )+sum(x(m+l  :m+N/2)  )  )  /N; 
xf=x; 


A. 7  S curve  Funtion  File 

S curve  calculates  the  S-Curve  for  an  input  sequence  as  defined  in  Chapter  III.  S curve  cor¬ 
relates  an  input  signal  with  a  locally  generated  signal,  which  is  also  an  input  to  the  function. 

This  function  allows  for  variable  correlator  spacing  so  that  narrow  or  standard  correlator  NCDLL 
S-curves  can  be  calculated  using  the  same  function  file.  The  Matlab  function  file  follows: 

function  [sc ,t]=scurve(x,c ,f s, s) 

7«  SCURVE  computes  the  tracking  curve  (s-curve)  for  a  non-coherent  delay  lock  loop  where  x 
7o  is  the  input  signal,  c  is  the  correlation  signal,  fs  is  the  sampling  frequency, 

7o  and  s  is  the  correlator  spacing.  If  s  is  not  specified,  the  default  is  1  chip. 

7  If  two  output  variables  are  used,  a  time  scale  for  the  s-curve  is  also  returned. 

7 

7  [sc ,t]=scurve(x,c ,f s , s) 

7 

7  Written  by  Chuck  Ormsby 
7  August  19,  1997 

7 

7  This  function  calls  the  vector_shift  subroutine: 
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% 


7o  Note:  Signal  names  used  in  this  file  refer  to  the  names  in  fig®  4-9,  pg  165  of 
7  Peterson,  Ziemer,  and  Borth:  Introduction  to  Spread  Spectrum  C<Mnications 
if  nargin==3 
s=l; 

end 

if  rem(s*fs,2)==0 
space=s*f s/2; 
else 

space=(s*f s+l)/2; 

end 

7o  Test  the  length  of  the  correlation  signals 
len=length(x) ; 
if  len>=50000 

trun=5e4;  %  For  very  long  sequences,  the  input  signal  must  be  tunc  ate  d  to  decrease 
7o  run  times 

else 

trun=length(x) ; 

end 

end 

x=x(l :trun) ; 
c=c(l :trun) ; 

7o  Produce  the  early  and  late  codes 
early=vector„shift (c , space) ; 
late=vector_shift(c , -space) ; 

7a  Produce  the  early  and  late  correlator  outputs 
yl=xcorr(early,x) ; 
y2=xcorr(late ,x) ; 

7o  Square  the  correlator  outputs 

zl=yl.^2; 

z2=y2.^2; 

7a  Low  pass  filter 
wcut=0.035; 

[b , a] =but t  er ( 5 , wcut ) ; 
zllp=f ilter(b,a,zl) ; 
z21p=f ilter(b,a,z2) ; 

7e  Create  the  s- curve 
sc=zllp-z21p; 

sc=sc  ./inax(sc) ;  7o  Normalize  the  s-curve 
t=linspace(-length(x) ,length(x) ,length(sc) )/f s ; 


A. 8  Vector-Shift  Function  File 


The  function  file  vectorshift  produces  the  time  shift  necessary  for  the  e^  and  late  locally 
generated  codes  in  the  NCDLL,  This  function  assumes  the  input  vector  is  one^iod  of  a  periodic 


92 


spreading  code.  The  function  then  shifts  the  elements  of  the  vector  circularly  left  or  right  to  create 
an  advanced  or  delayed  version  of  the  input  vector.  The  Matlab  function  file  follows: 
function  yp=vector_shift2(y,s) 

t  VECTOR.SHIFT  shifts  the  entries  in  the  vector  y  by  some  advance  or  delay.  A  delay  is 
%  entered  as  a  negative  number  while  an  advance  is  entered  as  a  positive  number. 

I 

7o  yp=vector_shift(y,s) 

I 

%  Note:  The  vector  y  is  assumed  to  be  periodic,  with  the  elements  of  y  comprising  an  integer 
7  number  of  periods. 

7o 

s=mod(s , length (y) ) ; 
if  s>=0 

right=y(l:abs(s)) ; 
left=y(abs(s)+l :length(y)) ; 
else 

right=y(l:length(y)+s) ; 
left=y(length(y)+s+l :length(y)) ; 

end 

yp=[left  right]; 


A. 9  Trackpoint  Function  File 

The  trackpoint  function  file  searches  for  the  tracking  point  in  a  vector  representing  the  NCDLL 
S-curve.  This  function  searches  for  the  point  where  the  S-curve  crosses  the  S{S)  —  0  axis  between 
the  minimum  and  maximum  points  of  the  S-curve.  The  function  then  reports  the  index  of  the 
data  element  of  the  zero  crossing  which  can  be  related  to  a  time  delay  through  knowledge  of  the 
sampling  frequency.  The  Matlab  function  file  follows: 
function  z=trackpoint (x) 

78  Trackpoint  searches  for  the  tracking  point  of  an  input  S-curve  vector.  The 
78  tracking  point  is  the  point  where  the  S-curve  crosses  the  S=0  axis  between 
78  the  minimum  and  maximum  values  of  the  S-curve.  If  the  vector  does  not  contain 
78  an  element  which  is  identically  zero,  a  linear  interpolation  between  the  two 
y,  nearest  elements  is  performed. 

78 

78  z  =  trackpoint  (x) 

78 

78  Written  by  Chuck  Qrmsby 

I 

78  This  function  calls  the  subroutines: 
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7o  maxfind.m,  minf ind . m ,  and  interpolate .m 

% 

inax_index=maxf  ind(x)  ; 
min_index=minf ind(x) ; 
if  min_index>max_index 
St art=max_ index ; 
stop=rain_index ; 
count=l ; 

for  index=start : stop 

if  x(index-l)>=0  &  x(index)<=0 

a( count , : )= [index- 1  index] ; 
count=count+l ; 

end 

end 


else 

start=min_ index ; 

St op=max_ index ; 

count=l ; 

for  index=start :stop 

if  x(index-l)<=0  &  x(index)>=0 
a( count , : )= [index- 1  index] ; 
count=count+l ; 
end 


end 

end 

[m,b]=interpolate(a(l) ,x(a(l)) ,a(2) ,x(a(2))) ; 
z=-b/in; 


A.  10  Maxfind  and  Minfind  Function  Files 

The  maxfind  and  minfind  function  files  search  for  the  mcLximum  and  minimum  points,  re¬ 
spectively,  in  an  input  vector.  These  functions  report  the  index  of  the  maximum  and  minimum 
values  of  the  input  vectors.  The  Matlab  function  files  follow: 

function  index=maxf ind(x) 

7  MAXFIND  returns  the  index  associated  with  the  maximum  value  of  a  vector. 

7o 

7#  index=maxf  ind(x) 

7c 

7c  Written  by  Chuck  Ormsby 
7c  July  24,  1997 
7o 

maxi=max(x) ; 
for  count=l : length (x) 
if  x(count)==maxi ; 
index=count ; 
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end 


end 


function  index=minf ind(x) 

%  MINFIND  returns  the  index  associated  with  the  minimuin  value  of  a  vector. 

I 

7o  index=ininf ind(x) 

7c 

7c  Written  by  Chuck  Ormsby 
7c  October  1,  1997 
7c 

mini=inin(x) ; 
for  count=l : length (x) 
if  X ( count )==inini; 
index= count; 

end 

end 


A.  11  Interpolate  Function  File 

The  interpolate  function  file  is  used  in  conjunction  with  the  trackpoint  function  to  determine 
the  tracking  point  of  an  input  S-curve  vector.  Interpolate  returns  the  slope  and  y-intercept  of  a 
line  joining  two  input  data  points.  This  file  is  used  by  the  trackpoint  function  to  calculate  the  slope 
and  y-intercept  of  a  line  joining  two  adjacent  data  points,  one  negative  and  one  positive,  in  the 
S-curve  vector.  From  knowledge  of  the  slope  and  y-intercept  of  this  line,  the  exact  zero  crossing 
point  can  be  calculated.  This  linear  interpolation  is  justified  because,  as  shown  in  Chapter  III,  for 
small  tracking  errors,  the  S-curve  is  approximately  linear.  The  Matlab  function  file  follows: 

function  [m, b] =interpolat e (xl , yl , x2 , y2 ) 

7c  INTERPOLATE  returns  the  slope  and  y-intercept  for  a  line  connecting  the  points  (xl,yl)  and 
7o  (x2,y2). 

7c 

7c  [m  b]=interpolate(xl,yl,x2,y2) 

7c 

m=(y2-yl)/(x2-xl); 

b=yl-m*xl; 
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A.  12  Summary 


This  Appendix  presents  the  Matlab  function  files  used  for  simulation  of  multipath  mitigation 
using  cepstral  techniques.  All  function  files  used  for  this  thesis  are  included  in  this  Appendix.  Any 
functions  not  specifically  included  are  built-in  Matlab  (version  5)  functions.  The  intent  of  including 
the  function  files  is  to  give  the  reader  a  better  understanding  of  the  function  of  the  cepstral  process 
for  multipath  mitigation.  Additionally,  by  including  these  function  files,  a  reader  familiar  with 
Matlab  should  be  able  to  reproduce  the  work  of  this  thesis. 
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